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Abstract 

The  accurate  and  computationally  efficient  estimation  of  signals  in  noise  has  long  been  a  field 
of  intense  study.  The  signal  present  in  natural  processes  is  many  times  well  modeled  as  the  sum 
of  real  or  complex  exponential  functions.  The  noise  for  computational  simplicity  is  often  assumed 
to  be  white  or  uncorrelated.  There  exist,  however,  many  cases  where  noise  is,  in  fact,  correlated. 
Accurate  and  efficient  estimates  of  the  signal  in  these  cases  require  that  the  noise  correlation  be 
taken  into  acconnt.  This  is  case  for  the  specific  application  of  interest  in  this  dissertation,  Synthetic 
Aperture  Radar  (SAR),  whose  images  of  objects  may  be  modeled  as  the  sum  of  two-dimensional 
complex  exponentials  (the  electromagnetic  scattering  centers  on  the  target). 

The  maximum  likelihood  estimate  of  the  signal  is  often  considered  the  best  possible  estimate  of 
the  signal.  While  many  white  and  colored  noise  maximum  likelihood  estimates  have  been  developed, 
efficient  solutions  to  the  estimation  of  one-  and  two-dimensional  exponentials  in  unknown  colored 
noise  do  not  exist. 

This  dissertation  develops  techniques  for  estimating  exponential  signals  in  unknown  colored 
noise.  The  Maximum  Likelihood  (ML)  estimators  of  the  exponential  parameters  are  developed. 
Techniques  are  developed  for  one  and  two-dimensional  exponentials,  for  both  the  deterministic 
and  stochastic  ML  model.  The  techniques  are  applied  to  Synthetic  Aperture  Radar  (SAR)  data 
whose  point  scatterers  are  modeled  as  damped  exponentials.  These  estimated  scatterer  locations 
(exponentials  frequencies)  are  potential  features  for  model-based  target  recognition. 

The  estimators  developed  in  this  dissertation  may  be  applied  with  any  parametrically  modeled 
noise  having  a  zero  mean  and  a  consistent  estimator  of  the  noise  covariance  matrix.  ML  techniques 
are  developed  for  a  single  instance  of  data  in  colored  noise  which  is  modeled  in  one  dimension  as  1) 
stationary  noise,  2)  autoregressive  (AR)  noise  and  3)  autoregressive  moving-average  (ARMA)  noise 
and  in  two  dimensions  as  1)  stationary  noise,  and  2)  white  noise  driving  an  exponential  filter.  The 
classical  ML  approach  is  used  to  solve  for  parameters  which  can  be  decoupled  from  the  estimation 
problem.  The  remaining  nonlinear  optimization  to  find  the  exponential  frequencies  is  then  solved 
by  extending  white  noise  ML  techniques  to  colored  noise.  In  the  case  of  deterministic  ML,  the  com¬ 
putationally  efficient,  one  and  two-dimensional  Iterative  Quadratic  Maximum  Likelihood  (IQML) 
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methods  are  extended  to  colored  noise.  In  the  case  of  stochastic  ML,  the  one  and  two-dimensional 
Method  of  Direction  Estimation  (MODE)  techniques  are  extended  to  colored  noise.  Simulations 
show  that  the  techniques  perform  close  to  the  Cramer-Rao  bound  when  the  model  matches  the 
observed  noise. 

Application  to  SAR  data  first  requires  that  damped  exponentials  have  not  been  distorted  by 
SAR  processing.  Then,  1-D  colored  noise  techniques  provide  better  estimates  at  low  model  orders 
(number  of  exponentials)  than  white  noise  techniques.  The  2-D  techniques  based  on  the  colored 
noise  model  also  more  accurately  model  SAR  data  than  existing  2-D  white  noise  techniques.  With 
an  appropriate  focusing  technique  and  matching  technique  for  the  exponentials  in  each  dimension, 
scatterers  are  located  with  high  resolution  in  SAR  images  and  colored  noise  techniques  improve 
these  location  estimates. 


Maximum  Likelihood  Estimation  of  Exponentials  in  Unknown  Colored  Noise 
for  Target  Identification  in  Synthetic  Aperture  Radar  Images 

I.  Introduction 

1.1  Background 

This  dissertation  presents  several  new  techniques  for  the  estimation  of  sinusoidal  or  exponen¬ 
tial  signals  in  the  presence  of  unknown  colored  noise.  They  include  techniques  for  estimation  in  one 
and  two  dimensions.  These  new  techniques  can  be  used  to  estimate  the  locations  of  electromagnetic 
scattering  centers  on  objects  viewed  with  Synthetic  Aperture  Radar  (SAR).  They  are  quite  general 
and  can  be  applied  to  any  exponential  estimation  problem  involving  colored  noise.  The  use  of  these 
estimation  techniques  although  centered  on  the  particular  characteristics  of  SAR  (a  single  instance 
of  data,  modeled  as  the  sum  of  damped  exponentials),  are  not  limited  by  these  characteristics,  but 
can  be  easily  extended  to  multiple  data  instances  and  undamped  exponentials. 

The  use  of  image  data  to  detect  and  identify  targets  is  pervasive  in  Air  Force  operations.  One 
particular  sensor  of  importance  due  its  stand-off  and  resolution  capabilities  is  Synthetic  Aperture 
Radar  (SAR).  A  SAR  image  is  formed  from  radar  target  returns  collected  over  a  flight  profile 
that  simulates  a  very  large  aperture  antenna.  This  SAR  image  data  is  used  by  aircrew  or  ground 
observers  to  provide  reconnaissance,  targeting,  or  navigation  information.  Attempts  to  convert  the 
SAR  data  to  a  reduced  set  of  derived  features  are  ongoing  efforts.  With  this  reduced  set  of  data 
more  sophisticated  pattern  recognition  algorithms  can  be  used  to  efficiently  detect  and  identify 
targets.  The  processor  required  to  create  SAR  images  with  sufficient  resolution  is  a  major  cost  in 
SAR  systems.  The  image  resolution  and  size  required  for  target  recognition  can  easily  exceed  the 
cost  associated  with  human  interpretation  of  the  data.  One  means  of  avoiding  the  cost  associated 
with  additional  resolution  is  to  directly  use  the  radar  scattering  data  or  complex  phase  history 
without  generating  a  SAR  image.  Pattern  recognition  algorithms  can  then  be  applied  directly  to 
features  extracted  from  the  phase  history  thereby  avoiding  the  cost  of  SAR  imagery. 

The  most  significant  of  these  extracted  features  provide  a  direct  correlation  to  the  locations 
of  electromagnetic  scattering  centers  on  the  target.  These  locations  can  be  estimated  with  great 
accuracy  by  parametrically  modeling  the  SAR  data.  To  efficiently  estimate  these  SAR  model 
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parameters,  the  model  must  conform  to  or  fit  the  SAR  data.  Methods  of  fitting  damped  exponentials 
to  SAR  data  can  more  accurately  estimate  scatterer  locations  when  unknown  background  clutter 
or  colored  noise  is  also  considered.  The  colored  noise  techniques  developed  in  this  research  improve 
the  accuracy  in  estimating  exponential  scatterer  locations  and  thus  increase  the  suitability  of  using 
exponential  scatterer  locations  as  pattern  recognition  features. 

1.2  Dissertation  Outline 

This  dissertation  is  developed  along  the  following  lines.  Chapter  2  introduces  the  problem 
solved  by  this  dissertation,  explores  the  previous  work  in  this  area,  and  describes  the  approach  in 
this  dissertation  to  solve  the  problem.  Chapter  3  presents  notation  and  mathematical  preliminaries 
used  as  a  basis  for  subsequent  developments.  Chapter  4  introduces  the  specific  problem  of  para¬ 
metric  estimation  of  scatterers  in  SAR  data.  The  assumptions  that  lead  to  the  damped  exponential 
model  are  examined.  Some  of  the  potential  pitfalls  that  are  associated  with  collecting,  processing, 
and  imaging  data  containing  damped  exponentials  are  presented.  In  Chapter  5,  the  parametric 
estimation  of  one-dimensional  damped  exponentials  in  colored  noise  is  examined.  Because  of  its 
robust  statistical  properties,  the  Maximum  Likelihood  (ML)  estimation  technique  is  applied  to 
solve  the  estimation  problem.  This  chapter  also  develops  several  new  estimation  techniques  that 
model  the  colored  noise  as  a  stationary  noise  sequence,  an  autoregressive  (AR)  noise  sequence,  or  an 
autoregressive  moving  average  noise  sequence  (ARMA).  Both  the  deterministic  ML  model,  which 
models  the  exponential  amplitude,  frequency  and  phase  as  unknown  constants,  and  the  stochastic 
ML  model,  which  models  the  amplitudes  as  random,  are  examined.  Chapter  6  examines  the  expo¬ 
nential  estimation  problem  in  colored  noise  for  two  dimensions.  The  1-D  deterministic  techniques 
developed  in  the  previous  chapter  are  extended  to  2-D  for  the  deterministic  ML  model.  Compu¬ 
tationally  less  intensive  techniques  that  estimate  each  dimension  independently  are  developed  for 
the  stochastic  ML  model.  Chapter  7  comes  full  circle  and  applies  the  colored  noise  techniques 
developed  to  1-D  and  2-P  SAR  data.  The  efficiency  of  the  colored  noise  techniques  in  fitting  the 
model  to  the  data  is  compared  with  the  white  noise  ML  techniques  and  other  techniques  such  as 
overmodeling  (fitting  more  exponentials  than  expected  to  the  data  and  selecting  those  that  best  fit 
the  data).  The  suitability  of  scatterer  locations  as  a  pattern  recognition  feature  is  also  examined. 
Chapter  8  presents  conclusions  and  some  potential  areas  of  future  research. 
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1.3  Contributions 


The  contributions  in  this  dissertation  include  the  solution  of  the  1-D  maximum  likelihood 
problem  of  estimating  deterministic  exponentials  in  unknown  colored  noise.  Computationally  effi¬ 
cient  estimation  techniqnes  for  several  noise  models  are  developed.  Similar  techniques  are  developed 
for  the  stochastic  maximum  likelihood  problem.  A  new  spectral  model  for  2-D  noise  is  developed 
leading  to  new  2-D  estimation  techniques  for  deterministic  and  stochastic  maximum  likelihood. 
Additionally,  a  mathematical  model  for  representing  damped  exponentials  on  an  irregularly  sam¬ 
pled  grid  is  developed  and  leads  to  a  computationally  efficient  method  for  interpolating  or  focusing 
SAR  images  containing  damped  exponentials. 
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11.  Problem  and  Approach 


2.1  Problem 

The  general  problem  solved  in  this  dissertation  is  most  easily  understood  in  terms  of  the 
array  processing  problem.  In  the  array  processing  problem,  the  aspects  of  exponential  estimation 
in  several  dimensions  can  be  observed.  In  addition,  because  of  the  general  nature  of  the  models  used 
in  this  dissertation,  the  approaches  derived  have  applications  in  all  areas  where  spectral  estimation, 
system  identification,  or  harmonic  retrieval  play  a  role.  The  parameters  of  any  process,  that  can 
be  modeled  as  the  sum  of  exponential  signals  in  unknown  Gaussian  noise,  may  be  estimated  with 
the  techniques  derived  in  this  dissertation. 

In  the  array  processing  problem,  a  set  of  sensors  are  placed  at  fixed  locations  in  a  wavefield 
that  consists  of  a  small  number  of  waves  as  shown  in  Figure  1.  These  sensors  then  measure  the 
amplitude  and  phase  of  the  wavefield  at  those  locations.  When  the  propagation  frequency  of  the 
waves  is  known,  an  array  may  be  designed  such  that  measurement  data  can  be  used  to  estimate  the 
directions  of  propagation  of  the  waves.  These  directions  of  arrival  (DOAs)  are  one  of  the  parameters 
estimated  by  the  techniques  of  this  dissertation. 

2.1.1  Wavefield  Model.  The  wavefield  model  whose  parameters  are  to  be  estimated  is 
now  constructed  for  the  specific  case  of  interest  in  this  dissertation,  the  electromagnetic  wavefield 
in  free  space.  Applications  to  many  other  wavefields  and  medium  exist,  including  for  example, 
acoustic  waves  in  the  ocean  (sonar)  and  acoustic  waves  in  rock  (seismography).  The  array  processing 
problem  for  the  electromagnetic  wavefield  is  quickly  derived  from  Maxwell’s  equations  [14]  [32].  The 
parametric  model  of  the  signal  measured  at  a  sensor,  at  time  t,  and  location  r  =  [  r*  Vy 
is 

y(r,t)  =  (1) 

where  s(r,t)  is  the  amplitude  of  the  wave,  w  is  the  propagation  frequency  of  the  wave,  and  k  = 
7[  cos  0  sin  sin  0  sin  0  coscj)  ]^  where  c  is  the  speed  of  light.  Now,  with  a  set  of  appropriately 
placed  sensors  the  wavefield  can  be  measured,  and  the  parameters  estimated  from  the 

measured  data.  The  position  of  the  sensors  samples  the  wavefield  at  various  locations,  for 
m  =  0..M  —  1,  and  the  sensors  are  seimpled  with  sampling  period  T,  at  discrete  times,  t  =  nT 
for  n  =  0..N  —  1,  to  produce  the  measured  data,  y(m,n).  The  simplest  array  geometry  (sensor 
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Figure  1.  A  Uniform  Rectangular  Array 


positioning)  is  given  by  the  uniform  linear  array  (ULA).  In  the  ULA,  the  sensors  are  placed  in  a 
line  at  a  fixed  interval.  With  the  uniform  linear  array,  only  one  direction  of  arrival  angle,  9,  can 
be  estimated.  A  simple  extension  of  the  ULA  to  two  dimensions  is  the  uniform  rectangular  array 
(URA),  a  rectangular  grid  of  sensors  r^imj  for  mi  =  O..M1  -  1,  m2  =  O..M2  —  1,  which  allows  both 
direction  of  arrival  angles,  6  and  <p,  to  be  estimated.  The  array  shown  in  Figure  1  is  such  an  array. 


2.1.2  One- Dimensional  and  Two-Dimensional  Data  Characterizations.  Depending  now 
on  the  statistical  assumptions  made  about  the  parameters,  the  measured  data  may  be  character¬ 
ized  in  one  of  several  one-  and  two-dimensional  estimation  problems.  Three  of  the  most  common 
characterization  are  detailed  here.  If  sensor  data  measured  a  different  times  is  considered  uncor¬ 
related,  the  data  is  temporally  white  and  the  sensor  data  for  each  time  sample  is  considered  a 
different  statistical  instance  of  the  sensor  data.  In  the  1-D  ULA  this  leads  to  a  one-dimensional 
estimation  problem  (parameters  {s,9,  known  w})  with  the  sensors  for  m  =  0..M  —  1  where 
the  data  from  each  time  sample  is  used  to  collect  statistical  information  (mean,  covariance)  on  the 
measured  sensor  data.  With  the  same  assumption,  the  2-D  uniform  rectangular  can  be  used  to 
solve  the  two-dimensional  estimation  problem  (parameters  {s,9,<j),  known  w})  with  sensors  r^imj 
for  mi  =  O..M1  —  l,m2  =  O..M2  —  1,  and  data  instances  n  =  0..N  —  1.  This  2-D  array  processing 
problem  may  also  be  extended  to  synthetic  aperture  radar  (SAR)  as  shown  in  Chapter  4.  The  last 
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data  characterization  considered  in  this  dissertation  is  that  involving  a  single  sensor.  This  case 
reduces  to  a  one-dimensional  estimation  problem  (parameters  {s,  with  a  single  data  instance, 
m  =  1.  This  is  a  time-series  analysis  or  one-dimensional  spectral  estimation  problem. 

As  always  in  estimation  problems,  the  goal  of  the  estimation,  accuracy,  is  limited  by  the 
available  data.  In  many  problems  the  numbers  of  sensors  is  extremely  limited.  A  simple  sonar 
array  may  contain  less  then  ten  sensors.  The  time  available  to  estimate  the  parameters  is  also 
limited,  both  from  signal  duration  and  reaction  time  (i.e.  as  in  a  fire  control  system),  and  drives 
the  computationally  complexity  of  the  estimation  technique.  Many  fast  estimation  techniques  have 
limited  accuracy  and  make  simplifying  assumptions  about  the  data  for  computationally  tractable 
results.  The  problem  then,  succinctly  stated,  is  to  accurately  estimate  the  parameters  of  particular 
array  and  data  model  chosen  in  a  computationally  efficient  manner  with  as  few  assumptions  on  the 
data  as  possible.  The  extension  of  the  array  processing  problem  to  the  scatterer  location  problem 
in  SAR  data  then  allows  the  application  of  the  array  processing  solutions  to  SAR  data. 

S.2  Previous  Work 

This  section  details  previous  research  that  has  been  completed  in  the  areas  that  this  disserta¬ 
tion  draws  on  as  a  basis  for  the  new  results  herein.  The  areas  covered  include  spectral  estimation, 
maximum  likelihood  estimation,  two-dimensional  exponential  estimation,  array  processing  in  col¬ 
ored  noise,  and  scatterer  identification  in  synthetic  aperture  radar  images.  As  always,  such  a 
summary  can  not  hope  to  be  complete  given  the  large  body  of  work  completed  in  each  of  the  areas. 
Historical  summaries  of  some  of  these  areas  are  found  in  [35]  [49]  [66] . 

2.2.1  Spectral  Estimation.  The  theory  involved  in  the  majority  of  the  research  in  spectral 
estimation  can  be  traced  back  to  Fourier  (1807).  These  nonparametric  techniques  are  appropriately 
named  Fourier  analysis.  It  is,  however,  Bunson  and  Kirchhoff  who  are  attributed  with  the  idea 
that  the  spectrum  of  the  light  (signal)  emitted  by  a  substance  can  be  used  to  characterize  that 
substance  and  its  physical  properties  [13].  The  inherent  utility  of  spectral  estimation  is  the  physical 
significance  of  the  representation  of  a  signal  by  its  spectrum,  whether  the  spectrum  represents  the 
molecular  content  of  a  substance,  the  key  vibrations  in  an  elastic  body,  the  stability  and  response  of 
a  dynamic  system,  or  the  frequency  and  direction  of  propagating  electromagnetic  waves.  Modern 
interest  in  characterizing  the  spectrum  of  random  (statistical)  signals  begins  with  Wiener  [83] 
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who  developed  the  autocorrelation  function  for  such  signals  and  its  Fourier  transform  relation 
with  the  power  spectral  density,  the  magnitude  squared  of  the  signal  spectrum.  Despite  robust 
performance  in  most  situations  the  Fourier  techniques  contain  inherent  limitations.  The  resolution 
of  the  Fourier  spectrum  is  limited  to  a  set  of  quantized  frequencies  and  the  variance  of  the  estimated 
spectral  energy  at  these  frequencies  is  equal  to  the  square  of  the  estimate.  Although  a  similar 
methodology  for  estimating  exponential  signals  can  be  traced  to  Prony  [55] ,  the  introduction  of  a 
parametric  model  by  Yule  [87]  to  characterize  random  processes  heralded  the  introduction  a  high 
resolution,  consistent  estimators  of  the  power  spectrum  of  a  random  process.  The  basis  for  both 
these  techniques  is  the  homogeneous  difference  equation  e[n]  =  0, 

apy[n]  +  ap-i2/[n  -  1]  H - h  aoy[n  -p]  =  e[n]  (2) 

where  the  Uj  for  i  —  0..p  are  the  coefficients  of  the  difference  equation  and  ?/[n]  are  samples  of  the 
random  process.  The  solution  to  the  difference  equation  is  the  sum  of  p  damped  exponentials 

=  (3) 

«=i 

where  the  Si  are  the  exponential  amplitudes  and  the  A,-  =  are  the  damped  exponentials. 

With  these  parametric  estimation  techniques,  the  accuracy  of  the  estimates  is  limited  only  by  the 
noise  present  in  the  data.  The  key  difference  between  the  approaches  of  Prony  and  Yule  is  that 
Yule  envisions  the  residuals  remaining  after  applying  the  difference  equation,  e[n],  to  the  be  a  white 
noise  driving  force  that  is  filtered  by  the  difference  equation  to  produce  the  colored  noise  process, 
y[n],y{n—\],  ■  ■  •  j/[n  — p].  This  is  the  first  of  the  rational  polynomial  noise  models,  the  autoregressive 
(AR)  noise  model.  Prony,  conversely,  considers  the  exponentials  to  be  fixed  deterministic  signals 
and  the  residuals  an  independent  white  noise  process.  Both  methods  arrive  at  the  same  spectral 
estimate  thus  they  both  may  be  used  to  estimate  both  exponential  signals  and  colored  noise.  A  key 
theorem  that  relates  these  two  types  of  random  processes  is  due  to  Wold  [84]  and  states  that  any 
random  process  may  be  decomposed  into  two  uncorrelated  components,  a  completely  deterministic 
process  that  may  be  exactly  estimated  from  its  past  samples,  and  a  random  process  that  may  be 
modeled  as  white  noise  driving  a  linear  difference  equation. 
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2. 2. 1.1  Modern  Spectral  Estimation.  In  modern  spectral  estimation  the  data  is 
modeled  as  instances  of  a  random  vector,  y„  =  [  Q  i  .  •  •  yn,M-iY  ~  ~ 

and  the  estimation  of  the  parameters  of  the  modeled  data  is  accomplished  by  eigenanalysis  of  the 
covariance  matrix  (or  some  related  matrix)  of  the  data, 

Ryy  =  E{yy*},  (4) 

where  E{-}  is  the  expected  value  operation  and  (•)*  indicates  complex  conjugate  transpose.  In  the 
case  where  the  data  is  modeled  as  a  linear  combination  of  signals  and  noise, 

y  =  Gs  +  w,  (5) 

where  the  columns  of  the  matrix  G  characterize  the  signals,  s  is  a  vector  of  the  signal  amplitudes, 
and  w  is  zero-mean  Gaussian  noise.  The  model  for  the  covariance  matrix  of  the  data  is  then 

Ryy  —  GPG*  -h  Rvrw,  (6) 

where  P  =  E{ss*}  and  iZ^w  is  the  covariance  matrix  of  the  noise.  The  eigenanalysis  methods, 
which  begin  with  Pisarenko  [54],  were  proven  with  Schmidt’s  multiple  signal  classification  (MUSIC) 
method  [61]  which  characterize  the  spectrum  of  the  data  in  terms  of  signal  and  noise  eigenspaces. 
For  the  case  of  white  noise  where  iZww  =  <r^  I  and  is  the  noise  variance,  this  allows  the  noise 
portion  of  the  data  (the  M  —  p  smallest  eigenvalues)  to  be  easily  discarded.  Other  eigenanalysis 
methods  included  Estimation  of  Signal  Parameters  via  Rotational  Invariance  Techniques  (ESPRIT) 
[57]  which  measures  the  signal  with  one  set  of  sensors  then  measures  the  signal  with  a  second  identi¬ 
cal  but  shifted  (in  location)  set  of  sensors.  An  eigenanalysis  of  these  two  measurements  produces  an 
estimate  of  the  signal.  In  the  last  type  of  eigenanalysis  method,  the  principal  components  method 
[76],  a  form  of  Prony’s  method  is  used  and  the  number  of  signals  is  overestimated.  Thus,  as  noted 
earlier  the  estimate  models  both  signal  and  noise.  As  in  the  previous  methods,  an  eigenanalysis  of 
the  data  allows  the  noise  part  of  the  measurement  to  be  discarded.  In  the  next  section,  the  close 
relation  between  eigenanalysis  method  such  as  MUSIC  and  the  maximum  likelihood  estimation  is 
explored. 
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2.2.2  Maximum  Likelihood  Estimation.  Maximum  likelihood  estimation  is  based  on  find¬ 
ing  the  probability  density  function  (and  model  parameters)  that  best  fits  the  observed  data.  In 
estimation  problems,  maximum  likelihood  (ML)  estimation  produces  accurate  estimators.  Several 
notable  large-sample  properties  of  ML  estimation  have  been  proven.  ML  estimation  produces  con¬ 
sistent  estimates  (i.e.  as  the  number  of  data  samples  grows  the  value  of  the  estimator  approached 
the  true  value).  ML  estimation  has  the  best  possible  efficiency.  If  an  unbiased  estimator  exists  that 
attains  the  estimation  performance  bound,  the  Cramer- Rao  bound  (see  Appendices  C  and  D),  it  is 
the  ML  estimator.  The  ML  estimator  is  asymptotically  efficient.  As  the  number  of  data  samples 
grows  the  variance  of  the  estimate  approaches  the  Cramer- Rao  bound  (CRB).  Additionally,  Monte 
Carlo  simulations  have  shown  that  ML  estimators  perform  well  on  small  size  data  samples  as  well. 

2.2.3  Deterministic  and  Stochastic  Maximum  Likelihood.  The  method  of  maximum  like¬ 
lihood  was  originated  by  Fisher  (1920)  whose  Fisher  information  matrix  characterizes  the  CRB. 
Contributions  to  the  method  were  made  by  Cramer  (1946),  Rao  (1946),  and  Wald  (1944).  The 
application  of  maximum  likelihood  estimation  to  the  array  processing  problem  can  be  traced  to 
[3]  [4]  [26]  [80] .  Most  notably  maximum  likelihood  estimation  has  been  applied  with  two  distinct 
assumptions.  The  first  assumption  is  that  the  signals  are  nonrandom,  in  which  case  the  signal  is 
the  mean  of  the  data 

E{y}  =  Gs,  (7) 

and  the  variance  of  the  data  is  due  solely  to  additive  noise 

Ryy  =  Rww  (6) 

This  case  is  denoted  as  deterministic  maximum  likelihood.  The  second  assumption  is  that  the 
signal  (amplitudes)  is  random.  Then,  the  data  has  zero  mean  and  covariance 

Ryy  =  GPG*  -1-  Rww  (9) 

where  P  =  E{ss*}.  This  case  is  called  stochastic  maximum  likelihood.  In  a  classic  series  of  papers 
Stoica  and  Nehorai  [68],  [71],  [72]  demonstrate  that  for  a  fixed  number  of  sensors,  m  =  0..M  —  1, 
as  the  number  of  data  instances  grows,  N  —>■  oo,  stochastic  ML  produces  asymptotically  better 


estimates  than  deterministic  ML.  Also,  in  the  papers  the  MUSIC  spectral  estimation  is  shown  to 
be  asymptotically  equivalent  to  deterministic  ML  [68].  Alternatively,  in  the  case  where  N  =  I  and 
the  number  of  sensors  grows,  M  — >  oo,  (or,  equivalently,  M  =  1,  N  oo),  deterministic  ML  and 
stochastic  ML  asymptotically  produce  the  same  results  [79].  All  of  the  ML  techniques  cited  thus 
far  require  some  form  of  search  of  the  signal  space  to  find  the  columns  of  the  matrix  G  that  best 
fit  the  data.  A  much  more  computationally  attractive  approach  applicable  to  the  ULA  or  time 
series  problem  was  proposed  independently  by  Kumaresan,  Scharf  and  Shaw  [33] ,  and  Dressier  and 
Macovski  [6].  In  this  approach,  the  deterministic  ML  problem  is  solved  by  iterating  a  linearized 
solution  to  the  problem.  The  deterministic  ML  solution  in  this  case  is  a  quadratic  optimization. 
With  this  technique,  constraints  are  easily  imposed  to  attain  solutions  that  are  a  restricted  set 
of  the  damped  exponential  model  (i.e.  exponentials  on  unit  circle  in  complex  plane,  real  valued 
exponentials,  or  sine  waves),  and  increase  estimation  accuracy  [34]  [62].  This  technique  is  called 
Iterative  Quadratic  Maximum  Likelihood  (IQML).  The  eigenspace  equivalent  of  this  technique. 
Method  of  Direction  Estimation  (MODE),  was  developed  by  Stoica  and  Sharman  [69].  With  an 
appropriate  weighting  factor  MODE  achieves  the  same  accuracy  as  stochastic  ML  as  the  number 
of  data  instances,  grows  N  oo  [70]. 

2.2.4  Two-Dimensional  Exponential  Estimation.  Several  two-dimensional  exponential  es¬ 
timation  techniques  have  recently  been  developed  that  do  not  require  a  computationally  intensive 
search  of  the  2-D  frequency  plane.  These  techniques  either  perform  a  one- dimensional  estimation 
once  in  each  dimension  and  pair  the  results  (a  1-D  by  1-D  method)  or  they  simultaneously  estimate 
the  2-D  frequencies  (a  full  2-D  method).  The  first  technique  by  Rao  and  Kung  [2]  uses  an  eigenanal- 
ysis  of  the  measure  data  and  a  technique  similar  to  ESPRIT  to  attain  the  spatial  frequencies  of  the 
exponentials  in  each  dimension.  These  frequencies  are  then  paired  according  to  which  combinations 
best  match  the  amplitudes  of  the  measure  data.  The  second  technique,  Matrix  Enhanced  Matrix 
Pencil  (MEMP)  from  [25]  also  exploits  the  principle  of  ESPRIT  to  find  the  spatial  frequencies. 
This  method  also  requires  pairing,  however,  it  efficiently  estimates  the  2-D  exponential  frequencies 
by  forming  an  estimate  of  the  2-D  covariance  matrix.  A  technique  called  2D  Prony  [59]  is  an 
extension  of  the  1-D  Prony  technique.  This  method  uses  a  technique  of  overmodeling  (estimating 
more  frequencies  than  the  data  contains)  similar  to  the  principal  components  method  and  has  a  ro¬ 
bust  method  of  selecting  the  two-dimensional  exponential  frequency  combinations  that  produce  the 
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highest  energy  signal.  The  remaining  two  techniques  are  derived  from  the  two  maximum  likelihood 
approaches  to  the  problem,  stochastic  ML  and  deterministic  ML.  The  stochastic  ML  technique, 
2D  MODE,  was  developed  by  Li  and  Stoica  [38].  This  method  uses  the  1-D  MODE  technique 
and  treats  each  columns  or  row  as  a  separate  data  instance  when  estimating  the  frequencies  of  the 
columns  or  rows.  A  full  2-D  technique  was  developed  by  Clark  [11]  which  follows  the  approach  of 
the  1-D  IQML  method.  This  method  characterizes  the  linear  space  orthogonal  to  linear  space  of 
the  signal  and  also  parameterizes  the  problem  in  terms  of  a  known  noise  covariance  matrix.  Since 
the  general  maximum  likelihood  problem  formulation  is  easily  parameterized  with  the  covariance 
matrix  of  the  noise,  the  ML  techniques  are  also  the  most  readily  modified  for  the  case  of  unknown 
colored  noise. 

2.2.5  Array  Processing  in  Colored  Noise.  The  white  noise  assumption  in  the  array 
processing  problem  is  in  some  cases  invalid.  Although  consistent  estimates  are  still  attained  by  white 
noise  techniques  in  the  presence  of  colored  noise,  additional  estimation  accuracy  can  be  attained  by 
incorporating  the  noise  coloration  in  the  model  [78].  When  the  noise  covariance  matrix  is  known 
the  most  common  colored  noise  technique  is  called  whitening.  In  this  method  the  covariance  matrix 
of  the  data  is  modified  to  make  the  noise  component  uncorrelated,  then  a  white  noise  technique 
is  applied  to  the  data.  In  most  cases  the  effect  of  this  whitening  on  the  signal  is  not  taken  into 
account  and  the  signal  estimate  is  still  not  optimal.  For  the  ULA  or  URA,  the  known  colored 
noise  methods  developed  by  Clark  [9]  produce  efficient  estimates  for  the  one-  and  two-dimensional 
deterministic  ML  problem. 

For  the  case  of  unknown  colored  noise,  several  approximate  maximum  likelihood  methods 
have  been  developed.  Such  methods  developed  by  Le  Cadre  and  Wax  preform  a  search  over  the 
potential  signal  and  noise  eigenvalues  to  find  those  that  best  represent  the  data  covariance  matrix 
[36],  [81].  The  methods  of  Bohme  [5]  and  Friedlander  [17]  characterize  the  noise  covariance  matrix 
as  a  linear  combination  of  known  matrices,  then  search  for  the  ML  signal  and  noise  estimate. 
Other  unknown  colored  noise  techniques  include  the  instrumental  variable  approach  [68]  which 
estimates  the  spatial  noise  using  temporally  uncorrelated  data  samples,  and  an  iterative  Toeplitz 
covariance  matrix  method  [77]  which  estimates  the  ML  mean  (signal)  and  Toeplitz  covariance 
matrix  (stationary  noise).  The  previous  techniques  all  involve  collecting  multiple  instances  of  data 
to  form  the  data  covariance  matrix. 
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The  first  unknown  colored  noise  method  for  the  single  data  instance  was  developed  by  Kay 
and  Nagesha  [29],  [47],  [48].  This  method  is  an  approximation  of  the  ML  estimate,  however,  the 
single  data  instance  problem  is  of  special  interest  in  the  SAR  scatterer  location  problem.  Also, 
the  exact  solution  to  unknown  colored  noise  problem  for  stochastic  ML  has  been  developed  by  Ye 
and  Degroat  [86].  Although  it  also  involves  a  multidimensional  search  for  the  ML  estimates,  the 
solution  provides  insight  into  the  relationship  between  white  and  colored  noise  techniques  and  how 
computationally  efficient  white  noise  techniques  might  be  extended  to  colored  noise. 

2.2.6  Scatterer  Location  in  Synthetic  Aperture  Radar  Images.  The  application  of  ex¬ 
ponential  estimation  techniques  to  synthetic  aperture  radar  is  quite  recent.  The  first  application 
involves  enhancements  to  the  images  themselves.  When  a  small  set  of  scatterers  is  identified  in  a 
SAR  image,  a  simplified  and  more  appealing  SAR  image  can  be  produced.  To  this  end,  Gupta 
applies  a  form  of  Yule’s  method  in  each  dimension  [18],  and  Hua,  et  al.  apply  the  MEMP  method 
[24] .  These  estimation  techniques  are  seen  to  estimate  the  locations  of  the  dominant  scatterers  in 
the  SAR  images.  A  foundation  for  the  work  in  scatterer  location  was  developed  by  Sacchini  who 
develops  estimates  of  the  scatterer  locations  in  SAR  data  of  simple  objects  (e.g.  flat  plates)  [58]. 
Several  recent  results  have  be  published  by  Li  involving  the  application  of  Fourier  and  model-based 
methods  to  computer  generated  SAR  data.  In  these  results  the  Fourier  based  methods  are  the 
most  robust,  and  the  performance  of  model-based  methods  is  significantly  degraded  due  to  the 
significant  amount  non-exponential  signal  present  in  the  SAR  data  [39],  [40].  This  degradation  is 
also  noted  in  [53]. 

2.3  Approach 

As  noted  above,  data  from  the  intended  application,  SAR,  contains  energy  not  well  modeled 
as  exponential  signal.  This  dissertation  investigates  modeling  this  energy  as  unknown  colored  noise. 
The  SAR  scatterer  location  problem  and  the  array  processing  problem  are  related,  thus  unknown 
colored  noise  solutions  to  the  array  processing  problem  are  explored.  Of  the  methods  presented  that 
solve  the  array  processing  problem,  the  maximum  likelihood  methods  provide  the  most  accurate  and 
consistent  estimates  of  the  array  processing  model.  These  methods  allow  for  use  of  a  parameterized 
noise  covariance  matrix  and  thus  can  eliminate  the  limiting  white  noise  assumption  and  provide 
more  accurate  solutions.  The  approach  then  begins  with  these  methods. 
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Next,  of  the  array  processing  characterizations,  the  one-dimensional  time  series  (or,  equiv¬ 
alently,  for  the  ULA,the  single  data  instance  case)  has  been  researched  the  longest  and  is  the 
most  well  understood.  Many  signal  and  noise  models  exist  for  time  series  analysis  and  maximum 
likelihood  estimates  for  these  models  are  well  developed.  Since  the  intended  application,  SAR, 
involves  diverging  electromagnetic  waves,  the  signal  model  is  the  damped  exponential  model  and 
the  noise  models  are  taken  from  the  set  of  rational  polynomial  models.  As  noted  earlier,  the  deter¬ 
ministic  and  stochastic  ML  assumptions  do  produce  different  results,  thus,  unknown  colored  noise 
solutions  under  both  assumptions  are  explored.  The  results  are  compared  in  terms  of  the  intended 
application,  a  single  data  instance  of  damped  exponentials  in  colored  noise. 

The  intended  application  is  two-dimensional,  thus,  2-D  techniques  are  investigated.  The 
two-dimensional  deterministic  and  stochastic  ML  methods  also  allow  for  use  of  a  parameterized 
noise  covariance  matrix.  These  two-dimensional  techniques  are  also  well  related  to  their  one¬ 
dimensional  counterparts,  thus,  the  methodology  employed  in  the  1-D  colored  noise  techniques  is 
easily  extended  to  two  dimensions.  Additionally,  the  different  two-dimensional  assumptions  of  2-D 
IQML  (2-D  exponentials  at  distinct  frequencies  in  each  dimension)  and  2-D  MODE  (exponentials 
at  any  intersection  of  a  1-D  by  1-D  frequency  grid)  allow  development  of  techniques  with  radically 
different  computational  and  scatterer  location  performance. 

The  array  processing  problem  and  the  SAR  scatterer  location  problem  are  related  through 
a  data  interpolation  or  focusing  step.  The  unknown  colored  noise  array  processing  techniques 
may  be  applied  to  SAR  data  after  the  SAR  data  is  interpolated  to  a  uniform  rectangular  grid. 
Improvements  in  scatterer  location  and  the  fit  of  colored  noise  models  to  the  data  are  investigated. 
The  effects  of  the  focusing  method,  colored  noise  model,  the  deterministic  or  stochastic  assumption, 
the  1-D  by  1-D  or  full  2-D  exponential  assumption,  and  the  1-D  by  1-D  matching  technique  used 
are  also  investigated. 
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III.  Mathematical  Preliminaries 


This  chapter  presents  some  mathematics  that  provide  a  foundation  for  the  developments  of 
this  dissertation.  Many  references  exist  for  the  following  material  on  matrix  algebra  and  random 
vectors.  The  principal  references  used  in  developing  these  sections  are  [27]  [41]  [51]  [66]. 

3.1  Matrix  Algebra  in  Signal  Processing 

In  signal  processing  as  in  many  areas  of  engineering,  matrices  and  linear  algebra  provide 
an  important  tool  for  solving  problems.  This  section  introduces  some  of  mathematical  concepts 
that  are  exploited  in  solving  the  problems  explored  in  this  dissertation.  Let  A  and  B  be  complex 
matrices  then  the  following  notation  defines  specific  operations  on  these  matrices. 

is  the  transpose  of  A. 

A*  is  the  complex  conjugate  transpose  of  A. 

vec(A)  is  the  vector  formed  by  stacking  the  columns  of  A. 

A®  B  is  the  Kronecker  product  of  A  and  B. 

AQ  B  is  the  Hadamard  (Schur)  product  of  A  and  B. 
is  the  pseudoinverse  of  A. 

Pa  is  the  projection  onto  the  column  space  of  A. 

3.1.1  Singular  Value  Decomposition.  One  facet  of  matrix  algebra  that  has  become  of 
increasing  importance  in  solving  complicated  problems  is  the  Singular  Value  Decomposition  (SVD). 
In  the  SVD  any  complex  valued  n  x  m  matrix,  H,  may  be  decomposed  as  the  product  of  three 
matrices 

H  =  USV*,  (10) 

where  (•)*  indicates  complex  conjugate  transpose.  Here,  the  nxn  matrix  U,  and  the  mxm  matrix 
V  are  unitary  matrices  {UU*  =  /  or  U~^  =  U*)  containing  the  left  and  right  singular  vectors, 
respectively.  5  is  an  n  X  m  diagonal  matrix  whose  unique  real-valued  elements  called  the  singular 
values  are  commonly  ordered  according  to  decreasing  magnitude.  The  matrices  U  and  V  are  not 
unique,  however,  they  provide  a  valuable  insight  into  the  range  (x  =  Rnxi,^G  Rmxi) 

and  null  space  {6  :  Hd  —  0)  of  the  matrix  H.  If  the  rank  of  H  (the  number  of  basis  elements  for  the 
range  of  H)  is  r,  then  S  contains  r  nonzero  real  singular  values  and  the  first  r  columns  of  U  form 
a  basis  for  the  range  of  H,  and  the  last  m  —  r  columns  of  V  form  a  basis  for  the  null  space  of  H. 
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Efficient  low  rank  (k  <  r)  approximations  of  H  can  also  be  formed  by  setting  to  zero  the  smallest 
r  —  k  singular  values  in  S. 

3.1.2  Cholesky  Decomposition.  A  special  case  of  the  SVD  occurs  when  H  is  Hermitian 
(H  —  H*).  Then  n  =  m  and  the  eigenspace  decomposition  of  H  is 

H  =  VDV*,  (11) 

where  V  contains  the  eigenvectors  of  H,  and  D  is  diagonal  with  real  eigenvalues.  When  D  con¬ 
tains  no  negative  elements,  H  is  positive  semidefinite,  and  the  singular  value  decomposition  and 
eigenspace  decomposition  of  H  are  the  same  with  U  =  V  and  S  =  D.  When  the  diagonal  of  D  also 
contains  no  zero  values,  H  is  positive  definite  and  the  Cholesky  decomposition  of  H  is 


H  =  VDiD'^V*  =  R*Q*QR=  R*R,  (12) 

where  D^V*  =  QR  is  the  QR  decomposition  of  D^V* ,  where  Q  is  unitary  and  the  n  x  n  matrix 
R  is  upper  triangular. 

3.1.3  Quadratic  Solutions.  Often  the  problem  to  be  solved  involves  an  optimization 
and  SVD  is  again  a  useful  tool.  The  recurrent  optimization  problem  in  this  dissertation  is  the 
minimization  of  a  quadratic  form 


ininy*A*Ay  =  inin||Ay||2 


(13) 


where  H-Hj  indicates  the  2-norm,  A  is  nxm,  and  y  =[  j/q  _x^  ]^.  When  the  solution  is  constrained 
such  that  ||y||2  =  1,  the  solution  is  minimizer  of  the  Rayliegh  quotient 


min 

y 


y*A*Ay 

y*y 


(14) 


and  the  minimizing  y,  ymin,  is  the  eigenvector  associated  with  the  minimum  eigenvalue  of  A*  A. 
When  the  constraint  is  j/o  =  1  and  A* A  is  positive  definite  then 


(AM)-ic 
“  c*(AM)-ic’ 


(15) 
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where  c  =[  l  V ■  When  yo  =  1,  and  A* A  is  rank  deficient  by  one,  then  restating  the  mini¬ 
mization  leads  to  the  least  squares  problem 


min 

X 


b  I  A 


1 

— X 


=  min||b— j4x||2 


(16) 


where  b  contains  the  first  column  of  A,  and  A  contains  the  last  N  —  1  columns  of  A.  The  solutions 
to  this  problem  are  the  Moore-Penrose  pseudoinverses 


Xmin 


=  i+b  =  I 


{A*A)-^A*h 

A*{AA*)-% 


m  <  n 
n  <  m 


(17) 


where  (•)+  indicates  pseudoinverse.  The  general  pseudoinverse  of  A  is  constructed  with  the  SVD 
and  includes  the  case  where  A  has  neither  full  column  or  row  rank.  The  SVD  pseudoinverse  is 


A+  =  VD+U* 


(18) 


where  D+  is  constructed  by  inverting  the  non-zero  elements  of  D. 


3.1.4  Projections.  The  pseudoinverse  then  leads  to  another  type  of  matrix  that  simplifies 
complex  solution  formula,  the  projection.  A  projection  may  be  constructed  for  the  matrix  A  as 


Pa=< 


AA+ 

A+A 


m  <  n 
n  <  m 


(19) 


The  eigenvalues  of  the  projection  matrix  are  zero  or  one,  and  the  eigenvectors  associated  with  the 
eigenvalue  one  form  a  basis  for  the  range  of  A.  Thus  for  m  <  n 


PaA  =  A. 


(20) 


Consequently,  the  projection  is  idempotent, 


=  f’A. 


(21) 


16 


Another  potential  way  to  simplify  an  optimization  is  to  perform  the  optimization  in  the  space 
orthogonal  to  the  original  problem  formulation.  For  example, 


maxyP^y  =  nhny(7  -  P^)y.  (22) 

When  a  matrix  representation  for  the  orthogonal  space,  G,  can  be  found  such  that  A*G  —  0  and 
the  matrix  [A  G  ]  is  full  rank,  the  projection  onto  the  orthogonal  space  can  be  constructed, 

Pg  =  G{G*G)-^G*,  (23) 

and  the  projections  form  a  resolution  of  the  identity  matrix, 

Pg  +  Pa  =  I.  (24) 


Thus, 


miny(7  -  PA)y  =  minyPey. 
y  y 


(25) 


3.1.5  Signal  Processing  Matrices.  Two  structured  matrices  play  a  special  role  in  signal 
processing  and  define  orthogonal  subspaces.  When  evenly  spaced  samples  of  a  periodic  function 
are  collected  and  placed  in  a  vector  y,  they  may  be  modeled  with  an  n  x  p  Vandermonde  matrix 
of  the  periodic  components,  G,  and  a  vector  of  the  component  amplitudes  s. 


y  = 


1 

Ai 


1 

Ap 


s  =Gs, 


(26) 


A 


n-1 

1 


\n  — 1 


The  matrix  that  defines  the  space  orthogonal  to  the  matrix  G  is  an  n  —  p  x  n  Toeplitz  matrix  (the 
elements  along  each  diagonal  are  identical).  This  matrix  A  contains  the  coefficients  of  a  polynomial 
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a(A)  ==  ttpA^  +  . .  .aiA^  +  ao  whose  roots  are  {Ai, . . Ap}, 


Thus  A*G  =  0. 


^0  ^1  •  •  •  dp  0 


0  dQ  dl  •  ■  ‘  dp 


(27) 


3.1.6  Kronecker  and  Hadamard  Products. 
Kronecker  product  which  is  defined  as 


Another  useful  matrix  algebra  tool  is  the 


A®  B  = 


anB  ai^B  ■■■  auB 
a^iB  a22B  •  •  •  a2iB 

o-kiB  Ok^B  ■■■  OkiB 


(28) 


where  A  is  a  ib  x  /  matrix,  S  is  a  m  x  n  matrix,  and  A  0  15  is  a  mk  x  nl  matrix  that  contains 
blocks  of  the  matrix  B  individually  multiplied  by  each  element  of  A.  With  the  Kronecker  product, 
operations  can  be  performed  simultaneously  on  the  rows  and  columns  of  a  matrix.  For  example,  in 
the  matrix  product  AXBA'  ,  A  and  B  are  separable  operations  that  act  on  the  rows  and  columns 
of  the  matrix  A,  respectively.  When  the  columns  of  the  matrix  AXBX  are  stacked  in  the  vec 
(vectorize)  operation  [41],  the  result  is  conveniently  given  in  terms  of  the  Kronecker  product  as 


vec(AA5^)  -{B®  A)vec(A). 


(29) 


The  Kronecker  product  has  several  useful  properties 


(A  0  B){C  ®D)  =  (AC)  0  {BD)  (30) 

(A  +  S)  0  (C  +  D)  =  (A  0  C)  +  (A  0  T>)  +  (5  0  C)  +  (5  0  D)  (31) 

tr(A  0  5)  =  tr(A)tr(5)  (32) 
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and 


(c  ®  Dyp  =  C°P  ®  D°P 


(33) 


where  op  includes  transpose,  conjugate  transpose,  pseudoinverse  and  inverse  when  the  individual 
inverses  exist.  Lastly,  the  Hadamard  (Schur)  product  of  two  matrices  of  the  same  size  is  defined  as 


aiifeii 

«125i2 

■  ■  ai/bii 

a2ihi 

0,22^22 

■  ■  Cl2ll>2l 

dklbkl 

0‘k2bk2  ■ 

■  ■  o-kl^kl 

(34) 


3.2  Random  Vectors 

If  the  elements  of  the  vector  x  are  samples  of  a  random  process,  then  the  distribution  of  the 
vector  can  be  described  with  a  probability  density  function  (PDF).  In  the  case  where  the  samples 
of  the  length  N,  complex  vector  x  are  complex  circular  symmetric  Gaussian,  the  PDF  is 


/(X) 


1 

7r^\<T^R\ 


where  for  the  expected  value  operation,  E{-}, 


(35) 


E{x}  =  p, 


(36) 


and 

E{{x  —  p)*  {x  —  p)}  =  R.  (37) 

Additionally,  for  complex  data  the  real  and  imaginary  parts  are  assumed  to  be  independent  with 

E{(x-p)^(x-p)}  =  0  (38) 

to  attain  a  valid  PDF.  The  matrix  R  is  the  covariance  matrix  of  the  vector  x  and  has  several 
interesting  properties.  The  matrix  is  Hermitian  symmetric  and  positive  definite.  The  covariance 
matrix  R  may  be  estimated  in  several  different  ways.  The  most  common  method  is  by  forming  the 
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outer  product  of  several  instances  of  the  vector, 


R  =  XX\ 


(39) 


where  X  contains  the  instances  of  x,  X  =  [  xi  X2  xx  ]•  This  method  provides  a  consistent 
estimate  of  R  as  K  —y  oo,  and  when  K  >  N,  the  covariance  matrix  is  full  rank  with  probability 
1.  The  estimated  R,  however  is  not  positive  definite.  In  the  case  where  x  is  stationary,  a  positive 
definite,  full  rank,  covariance  matrix  can  be  formed  as  a  Toeplitz  matrix  containing  the  estimated 
autocorrelation  sequence  of  x,  f[i]  for  i  =  —{N  —  1) . .  .N  —  1, 


R  = 


f[0]  1] 

f[l] 

r[N  -  1] 


•  •  •  f  [1  -  AT] 

f[-l] 
f[l]  f[0] 


(40) 


This  is  equivalent  to  letting 


a:  = 


a:[0]  a:[l] 

0 

0 


x[N  -  1]  0 


*[0] 

0 


x[l] 

a;[0] 


x[l] 


x[N  -  1]  0 

a:[Ar  -  1] 


(41) 


be  a  JV  X  (2N  —  1)  matrix  in  Equation  39.  Now,  since  R  is  positive  definite,  its  Cholesky  decom¬ 
position,  R  =  L*  L,  provides  an  efficient  method  for  creating  realizations  of  noise  whose  covariance 
is  R.  The  noise  sequence  x  whose  covariance  is  R  is  constructed  from  L  as 


x  =  L*6 


where  6  is  vector  whose  elements  are  a  unit  variance  white  noise  sequence. 


3.2.1  Autocorrelation  Sequence  and  Power  Spectral  Density.  The  autocorrelation  se¬ 
quence, 

.  N-k 

^  a;[n]a;[n-|- fc]  =  — AT -f  1 . .  .0  . . . IV  -  1,  (42) 


n=0 
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is  efficiently  computed  via  the  Fast  Fourier  Transform  (FFT)  as  the  Discrete  Fourier  Transform 
(DFT)  of  the  data.  When  FFT  is  zero-padded  to  length  2N  —  1,  it  interpolates  all  2N  —  1  auto¬ 
correlation  coefficients  required  to  form  the  matrix  R.  This  Fourier  transform  of  the  data  can  also 
be  viewed  a  Fourier  transform  of  samples  of  the  Power  Spectral  Density  (PSD)  of  the  data, 


P.xien 


n=0 


(43) 


In  this  case,  the  data  samples,  a:[n]  n  =  0 . . .  —  1,  are  coefficients  of  a  filter  satisfying  a 

linear  difference  equation.  In  the  most  general  form,  this  linear  difference  equation  is  a  function  of 
both  the  input  sequence  and  the  output  sequence.  The  sequence  then  satisfies  an  autoregressive 
moving  average,  ARMA(o,  5),  difference  equation. 


6o®[n]  =  — 6ia:[n— 1]  . . .  —  6(,x[n  —  o] 

-l-Co€[n]  +  cie[n  -  1]  . . .  -|-  Cge[n  -  q], 


(44) 


where  e[n]  is  the  input  sequence.  This  allows  the  PSD  to  be  characterized  as  the  rational  function 
of  two  polynomials.  The  ARMA  PSD  is 


|ELoCne-^'"‘- 

\En=obne-i-‘^ 


(45) 


A  special  case  of  this  PSD  occurs  when  the  single  coefficient  the  numerator  is  cq.  This  produces 
the  autoregressive  (AR)  PSD. 


3.S.2  Maximum  Likelihood.  Methods  of  efficiently  estimating  the  parameters  in  these 
models  are  of  interest.  If  an  efficient  estimator  exists,  the  Maximum  Likelihood  (ML)  technique 
is  guaranteed  to  find  this  estimator.  In  this  case,  the  ML  estimator  is  unbiased  and  attains  the 
Cramer-Rao  estimation  bound  (CRB)  (see  Appendices  C  and  D).  The  ML  estimator  has  also  been 
shown  to  be  asymptotically  efficient  and  asymptotically  normally  distributed  (for  a  large  number 
of  data  samples,  high  signal-to-noise  ratio,  or  large  number  of  data  instances). 
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In  the  maximum  likelihood  technique,  the  PDF  of  the  data  is  evaluated  at  the  observed  data. 
Then,  the  unknown  parameters  of  this  PDF  function  are  optimized  to  find  values  of  the  parameters 
which  maximize  the  function.  These  values  are  the  ML  estimates  of  the  parameters. 
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IV.  Modeling  of  Synthetic  Aperture  Radar  Data  with  Exponentials 


This  research  develops  exponential  estimation  techniques  that  can  be  applied  to  data  sets 
such  as  Synthetic  Aperture  Radar  (SAR)  data.  As  a  foundation,  the  relationship  between  the 
damped  exponential  model  and  the  SAR  process  is  established  here.  When  SAR  data  is  collected 
under  certain  constraints,  the  data  can  be  modeled  with  a  damped  exponential  model.  A  spatial 
(k  =[  kx  ky  ]^)  and  temporal  frequency  (w)  domain  analysis  of  SAR  data  collection  shows  that 
SAR  data  may  be  used  directly  in  exponential  estimation  algorithms  with  the  addition  of  an  inter¬ 
polation  step.  In  the  interpolation  step,  the  k— w  domain  data  or  phase  history  is  transformed  from 
the  sampling  grid  on  which  it  was  attained  (polar  for  chamber  data,  or  polar-linear  for  airborne 
data)  to  a  rectangular  grid.  This  interpolation  step  may  be  skipped  provided  the  k  —  w  domain 
data  are  collected  on  a  grid  that  is  very  nearly  rectangular.  The  error  induced  by  this  approxi¬ 
mation  is  examined  later.  Exponential  estimation  techniques  model  SAR  data  directly.  Thus,  any 
modification  of  the  data  will  degrade  the  estimates.  The  use  of  filter  windows,  zero-padding,  or 
inappropriate  means  of  fitting  the  k  —  w  domain  data  to  a  rectangular  grid  significantly  degrades 
the  performance  of  exponential  estimation  techniques.  Like  exponential  estimation,  the  interpola¬ 
tion  step  also  involves  approximating  from  measured  data.  Thus,  interpolation  limits  estimation 
accuracy.  A  good  interpolation  step  or  an  estimation  step  that  does  not  require  interpolation  are 
of  interest  for  accurately  locating  scattering  centers. 

4.1  SAR  Data  Collection 

SAR  data  collection  is  modeled  here  in  a  two-dimensional  perspective  as  shown  in  Figure  2. 
Here  airborne  SAR  data  is  collected  by  spotlighting  a  ground  target  location.  For  a  linear  flight 
path  the  collection  geometry  involves  a  line  (the  flight  path)  and  a  non-colinear  point  (the  target). 
The  data  collection  can  be  thought  of  as  occurring  entirely  in  the  plane  represented  by  the  point 
and  the  line  (the  slant  plane).  The  target  signature  is  collected  as  if  the  target  is  viewed  orthogonal 
to  this  plane.  The  following  simplifying  assumptions  are  made;  1)  the  aircraft  flight  path  is  linear, 
2)  the  radar  transmits  and  receives  the  target  echo  at  the  same  point  in  space,  and  3)  the  Doppler 
in  the  signal  is  negligible.  This  model  can  also  be  expanded  to  include  deviations  from  the  linear 
path  to  account  for  actual  aircraft  flight  path,  or  a  circular  path  (chamber  data). 
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Synthetic 

Aperture 


Figure  2.  Spotlight  SAR  image  formation  geometry 

With  the  above  assumptions,  and  following  the  k  —  w  domain  analysis  of  [8]  [63]  [64]  [65] 
the  signal  received  by  the  radar  at  time  t  and  position  u,  s(u,t),  is  the  integral  of  the  transmitted 
pulse,  p,  and  the  slant-plane  target  reflectivity  function,  f{x,  y),  where  x  and  y  are  coordinates  in 
the  slant  plane.  Then 

s{u,t)=:  j  J  j  f{x,  y,  z)p{t  -  (a:  -  Xaf  +  {y-u-  yaY)dxdydz,  (46) 

where  c  is  the  speed  of  light,  {xa,ya)  is  the  center  of  the  synthetic  aperture,  ^  is  the  distance  from 
the  slant  plane,  and  the  limits  of  integration  in  this  and  all  subsequent  indefinite  integrals  are  ±oo. 
If  we  take  the  Fourier  transform  with  respect  to  t,  defined  as. 


/OO 

s{u,t)exp{-jut)dt, 

-OO 


and  filter  s{u,t)  with  the  matched  filter  for  the  pulse  p,  the  signal  becomes 


(47) 


P*(w)5(u,w)  =  |P(w)|^y  J  J  f{x,  y,  z)  exp{-j2k\/  {x  -  x^f  +  {y  -  u  -  yaY)dxdydz,  (48) 


where  k  =  and  P*{u))  is  the  response  of  the  matched  filter. 
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Next,  assume  that  the  target  consists  of  point  scatterers  at  distinct  locations  in  space.  Then 


f{x,  y,z)  =  'Y^  f„6{x  -  Xn)S{y  -  yn)S{z  -  z„),  (49) 

n 

where  the  amplitude,  includes  target  reflectivity,  and  a  damping  factor  due  to  wave  divergence 
which  is  inversely  proportional  to  (a;„  —  Xa)^  +  («/«  —  u  —  yaY-  The  image  perspective  orthogonal 
to  the  SAR  collection  plane  includes  scatterers  at  all  depths,  z,  in  the  target,  thus,  the  dependence 
on  z  integrates  out,  giving  a  received  signal 


P*{u})S{u,u)  =  \P{w)\^  j  J  f{x,  y)  exp{-j2k\J  {x  -  Xaf  +  (y  -  u  -  yaY)dxdy,  (50) 


where  f{x,  y)  -  ~  Xn)Kv  “  Vn)- 

The  Fourier  transform  with  respect  to  u  is  given  by  a  plane  wave  decomposition  of  the 
spherical  wave,  exp  ^-j2k^ {x  -  Xa)^  +  (y  —  u  —  ya^'j  ,  which  gives  the  Fourier  transform  relation 
[64] 

Pu{exp(-j2k\/x^  +  (y-  w)2)}  =  exp(-j\/4k^^^x  -  jkuy)  (51) 


where  the  Fourier  transform  with  respect  to  u,  Pu,  is  defined  as. 


/OO 

s{u)exp{—jkuu)du. 

•CO 

Then,  the  —  w  spectrum  of  the  received  data  is 


(52) 


P*{uj)S{ku,u)  =  |P(w)|^y  J  f{x,y)exp{-jy/^je^^^{x  -  Xa)- jku{y -ya))dxdy  (53) 

=  \P{u)\^exp(jy/ik^  -  k^Xa+jkuya)  J  J  f{x,y)exp{-j\/4:k^  -  klx  -  jkuy)dxdy 
=  |P(w)|^exp(j\/4A;2  -  klxa  dr  jkuya)F{\/^k^  -  kl,ku), 


where  F  is  the  Fourier  transform  of  /.  If  we  normalize  the  received  signal  by  l/|P(w)p,  it  is 
simply  the  2-D  Fourier  transform  of  the  target  reflectivity  function  f{x,  y)  with  a  phase  shift  term 
to  account  for  the  target  angle  with  respect  to  the  aperture  (squint).  If  samples  of  P(-,  •)  were 
available  at  uniformly  spaced  intervals  of  ^4fc^  —  k^  and  ku,  f{x,y)  could  be  recovered  error- 
free  by  this  procedure  [64].  The  data  collection  geometry,  however,  does  not  permit  this  and 
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uniformly  spaced  samples  must  be  interpolated  from  the  available  data.  The  most  common  form  of 
interpolation  is  to  assume  that  a  point  scatterer  exists  at  each  sample  in  the  image,  f(x,  y).  This 
is  called  focused  SAR. 

Focused  SAR  for  chamber  data  [42]  involves  phase  shift  terms  to  align  the  samples  in  lines 
perpendicular  to  the  kx  and  ky  axes..  The  focused  SAR  formulation  is 

1  ^TT  pCO 

f{x,y)  =  —  /  /  F{k,  6)  exp{—j2kx  cos  0  —  j2ky  sin  9)  k  dkdO  (54) 

Stt  Jo 

^2^  poo 

=  —  /  /  F{2k,6)exp(j2kx{l  —  cos  6)  +  jy(9  —  2k  sin  9))  exp{—jx2k  —  jy9)2dkd6. 

Stt  Jo  J-oo 

Fourier  transforming  both  sides  of  Equation  54  shows  that  the  collected  data  F[k,  9)  can  be  simply 
modulated  by  exp(j2jfca;(l  —  cos  9)  +  jy{9  —  2fcsin^))  to  attain  the  rectangularly  sampled  F(kx,  ky). 
In  this  method,  the  contribution  of  every  measurement  is  added  coherently  with  respect  to  each 
pixel,  assuming  that  a  point  scatterer  exists  at  every  pixel.  This  focusing,  however,  degrades 
damped  exponential  signals  in  the  data  by  convolving  them  with  a  point  spread  function,  and  limits 
the  ability  to  attain  high  resolution  estimates  of  the  exponentials.  More  appropriate  methods  of 
interpolation  are  discussed  in  the  following  sections. 

4.2  Fitting  the  Damped  Exponential  Model 

The  point  scatterer  assumption  in  the  SAR  process  leads  directly  to  the  exponential  model 
with  a  damping  factor,  /„,  to  account  for  wave  divergence.  With  the  point  scatterer  assumption, 
the  Fourier  transform  of  f{x,  y)  is  the  2-D  damped  exponential  model 

F{kx,ky)  =  Y^fnexp{-jkxX„)exp{-jkyyn).  (55) 

n 

Notice  that  this  it  equivalent  to  the  array  processing  model 

=  ^fn  exp(-iwt  -  jk^x„)  t  =  0,  (56) 

n 

where  k  =[  kx  ky  V  Xn=[  a;„  j/„  V  time  t  =  0.  Thus,  the  array  processing  model  for 

which  many  parameter  estimation  techniques  have  been  developed  can  be  exploited.  The  difSculty 
in  applying  this  model  lies  in  the  sampling  of  the  target  reflectivity  function. 
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airtome  sampling  grid  chamber  sampling  grid 


l<_x  (radians)  k_x  (radians) 


Figure  3.  Target  Spatial  Frequency  Samples 

The  calculation  of  the  Fourier  transform  of  the  target  reflectivity  function  (for  the  airborne 
case)  is  error  free.  However,  when  samples  are  taken  uniformly  in  u,  the  samples  are  not  uniform 
in  kx,  and  when  the  path  is  circular,  the  samples  are  on  a  polar  grid  in  the  {kx,ky)  plane.  The 
sampling  grid  for  the  airborne  and  chamber  cases  are  plotted  in  Figure  3.  Interpolation  must  be 
performed  between  the  data  sampling  grid  and  the  rectangular  grid  needed  to  estimate  damped 
exponential  parameters.  Several  focusing  methods  that  create  SAR  images  do  this.  However,  these 
assume  that  a  point  scatterer  exists  at  each  pixel  and  thus  impose  a  point  spread  function  on  the 
data.  This  windowing  degrades  the  high  resolution  estimates  of  parametric  estimation  methods 
[30] .  Inappropriate  means  of  interpolating  a  SAR  image  such  as  zero-padding  X- PATCH  data  also 
destroys  the  performance  of  exponential  estimation  techniques  [53].  Gradient  methods  also  exist 
for  focusing  SAR  data  at  existing  scatters,  but  these  techniques  are  computationally  burdensome 
[15]  [20].  This  leaves  three  solutions:  1)  collect  data  on  an  almost  rectangular  grid,  2)  interpolate 
the  data,  3)  estimate  exponentials  from  nonuniform  data. 

4.S.I  Rectangular  Assumption.  The  rectangular  assumption  has  been  used  in  SAR  tech¬ 
niques  that  incorporate  exponential  estimation  [24]  [58].  The  error  incurred  by  using  the  polar  or 
polar-rectangular  grid  data  directly  is  determined  by  examining  the  deviation  of  these  grids  from 
the  rectangular  grid  (see  Figure  4).  For  chamber  data  on  the  polar  grid,  the  data  are  collected 
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Normalized  Spatial  Frequency  Samples  -  Almost  Rectangular 
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Figure  4.  Normalized  Target  Spatial  Frequency  Samples  for  Chamber  SAR.  The  small  angular 
and  frequency  increments  cause  the  grid  to  be  almost  rectangular. 


about  an  angle,  <j>,  from  the  a;-axis,  in  a  beam  of  ^  €  [— ^ir])  and  at  frequencies  k  £  [feminj  ^max]- 
The  first  requirement  is  to  achieve  better  than  Fourier  resolution,  thus,  each  sample,  F{krrn,9m2), 
must  lie  in  a  separate  cell  in  the  rectangular  grid.  In  Figure  4  if  the  first  sample  in  a  row  is  lower 
than  the  last  sample  of  an  adjacent  row,  then  resolution  of  a  scatterer  to  less  than  a  pixel  (Fourier 
bin)  is  not  possible.  This  occurs  first  in  Figure  4  in  top  row  of  samples  when 


N -2, 

<  - 

N-1  ” 


(57) 


The  error  in  point  scatterer  location  is  given  by  the  frequency  estimation  accuracy.  This  accuracy  is 
bounded  by  the  maximum  kx  or  ky  deviation  of  samples  within  a  column  or  row  of  the  rectangular 
grid.  This  error  consists  of  two  orthogonal  components,  and  Cky.  Since  the  polar  data  are 
circularly  symmetric,  (j)  is  arbitrary  and  is  chosen  to  be  O'*.  The  error  in  the  k^  direction  is  greatest 
at  the  outer  ring  and  is 


^kx  —  ^max  COS  0  iTj^iax  COS  —  ^max(I  COS  0^). 


(58) 
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The  error  in  the  y  direction  is  the  deviation  of  the  radial  lines  from  parallel  beam  center  and  is 

City  =  ^max  sin  Ot  ^min  sin  Oh  =  (^max  ^min)  sin  Of,.  (^9) 

The  maximum  error  is  then  the  maximum  of  these  errors.  Airborne  SAR  is  regularly  sampled  in 
the  ky  direction  and  thus  contains  only  the  ejt^  error  term. 

4.S.S  Interpolation.  Since  the  k  —  w  domain  data  resembles  the  sum  of  complex  expo¬ 
nentials  the  sampled  data  is  quite  smooth.  Intermediate  samples  may  be  estimated  from  the  polar 
sampled  data  to  form  a  rectangular  grid.  The  simplest  form  of  interpolation,  phase-shifting  the 
data  to  a  rectangular  grid  was  described  earlier.  More  complex  methods  of  interpolation  including 
1-D  and  2-D  Lagrange  interpolation  polynomials  [37]  and  an  inverse  distance  method  [60]  are  also 
examined. 


4-2. 3  Scattering  Center  Focus.  This  section  introduces  a  new  compact  model  for  non- 
uniformly  sampled  2-D  exponentials.  With  this  model  a  new  computationally  efficient  method  of 
focusing  SAR  phase  history  data  by  simply  applying  phase  shifts  to  the  data  samples  is  developed. 

The  phase  history  data  F,  the  Fourier  transform  of  the  SAR  image,  can  be  written  as  2-D 
damped  exponential  data  from  p  point  scatterers  and  noise  N 

0  <  mi  <  Ml  —  1 

F[mi,m2]  =  +  A[mi,m2]  (60) 

i=i  0  <  m2  <  M2  —  1, 

where  the  Aj  and  ji  are  the  complex  spatial  frequencies,  and  the  Sj  are  the  complex  amplitudes  of 
the  scatterers  and  A[mi,  m2]  is  additive  noise.  Arranging  F  in  an  Mi  x  M2  matrix  gives 

F  =  GSH'^  -k  W,  (61) 

where 

^=[gl  g2  •••  gp  ]  gi  =  ^Mi(Ai), 

5  =  diag([si  S2  •••  Sp  ]), 

H  =  [hi  h2  •••  hp  ]  hi  =  $M,(Ti). 
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and 


^Af(^)  =  [  1  2:  ■■■  Y  ■ 

This  model  assumes  that  the  data  are  sampled  on  a  rectangular  grid.  A  reformulation  of  Equation 
61  will  expose  how  the  rectangular  grid  is  involved  in  the  model  and  how  it  may  be  replaced  with 
a  different  grid  to  reflect  the  data  collection.  Consider  a  2-D  damped  exponential  with  a  single 
exponential  in  each  dimension.  The  matrix  form  may  then  be  expressed  as 


Fz=SiXf&yy 


(62) 


where 


1 

1 

•  1 

1 

2  ■ 

•  M2 

2 

2  • 

•  2 

,  Y  = 

1 

2  • 

■  M2 

Ml 

Ml  • 

••  Ml 

1 

2  • 

•  M2  _ 

(63) 


and  the  scalar  Aj  taken  to  a  matrix  power  results  in  a  matrix  of  the  same  size  with  elements 
where  Xmirm  is  the  mim^^  element  of  X  and  ©  represents  the  element-wise  Hadamard  or  Schur 
product  of  two  matrices.  X  and  Y  represent  the  coordinates  of  the  sampling  grid  on  which  the  data 
is  taken.  To  express  a  damped  exponential  which  is  sampled  with  a  non-rectangular  grid,  simply 
replace  the  matrices  X  and  Y  with  the  rectangular  coordinates  of  that  grid.  For  the  multiple 
exponential  case, 

(64) 

i=l 


Each  element  of  F  is 


■  mifriQ 


*=i 

p 


i=l 


The  damped  exponential  model  is  distorted  at  each  pixel  by  the  weighting  factors 


ymimo  ^2'! 

\  Ti  /•  1  • 


For  conversion  to  a  normalized  polar  grid  these  weights  are 


m j  ' 


a,nd  0^2  =  m2^9 .  Assume  that  the  A,  and 
7i  are  undamped  exponentials.  Then,  m2  gjj^pjy  applies  a  phase  shift  to  each 
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term  in  Fmim^  average  phase  shift  for  equal  amplitude  scatterers  is 


-mim2 


,  COS  6 


^  t=l 


*»2  *”‘1  _|_  ' 


(65) 


Thus,  the  focusing  term  in  Equation  54  is  clearly  inappropriate  since  it  assumes  that  a  scattering 
center  exists  at  every  pixel  in  the  image.  This  focusing  distorts  scattering  center  locations. 

Equation  65  is  easily  implemented.  Estimate  the  Aj  and  7,  by  picking  the  brightest  pixels  in 
the  image  and  calculating  their  spatial  frequency.  Any  other  2-D  frequency  estimation  technique 
can  be  used.  However,  since  only  the  scattering  center  pixel  locations  are  required,  the  DFT  of 
F{kmi,Sm2)  provides  the  best  result.  This  method  also  provided  a  fast  means  of  improving  SAR 
image  quality  when  2-D  interpolation  methods  are  too  intensive.  Figures  5,  6,  and  7  show  images 
resulting  from  the  various  focusing  methods  used  on  chamber  data  of  a  C-29  aircraft  and  the 
unfocused  data  (rectangular  assumption). 


4.2.4  1-D  and  2-D  Interpolation.  For  the  1-D  Lagrange  interpolation  polynomial,  the  ky 
dimension  for  chamber  data  was  interpolated.  An  TV**  order  interpolation  polynomial  was  formed 
for  each  frequency  in  the  polar  grid.  This  polynomial  was  then  evaluated  at  the  rectangular  grid 
points  giving  the  interpolated  data.  The  most  computationally  intensive  techniques  examined  were 
2-D.  The  inverse  distance  method  used  the  2-D  gradient  to  interpolate,  and  a  Lagrange  polynomial 
interpolated  the  stacked  columns  of  the  2-D  data.  As  a  basic  measure  of  the  relative  merit  of 
each  technique,  an  ideal  image  of  point  sources  was  sampled  on  an  ‘almost  rectangular’  polar  grid, 
then  interpolated  by  each  technique  to  a  rectaugular  grid.  Figure  8  shows  the  energy  loss  of  each 
interpolation  method,  e  =  where  ||  •  ||f  indicates  the  Frobenius  norm.  As  expected  the 

2-D  interpolation  methods  show  less  degradation. 


4.3  Results 

Accuracy  is  important  in  determining  relative  scatterer  locations.  When  exponential  tech¬ 
niques  model  un-interpolated  data,  minimizing  the  energy  between  the  model  and  an  image  does 
not  accurately  locate  the  scatterers  in  an  image.  Thus,  ideal  data  is  used  here  to  determine  the  per¬ 
formance  of  interpolation  and  estimation  using  SAR  data.  Ideal  SAR  data  of  two  point  scatterers 
was  created  to  test  the  accuracy  of  several  exponential  estimation  techniques  with  different  interpo- 
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cross  range  (meters) 


Figure  5.  C-29  focused  with  the  technique  developed  by  Mensa  [42] 
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cross  range  (meters) 

Figure  6.  C-29  focused  with  scattering  center  focus  (64  scattering  centers) 
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Figure  7.  Unfocused  image  of  C-29  (rectangular  assumption) 
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Normalized  Error  in  Fitting  Rectangular  Grid 


Figure  8.  Interpolation  error  (Chamber  Data) 


lation  methods.  Mean  Square  Error  (MSE)  for  phase  history  data  generated  from  Equation  60  for 
two  exponentials  with  parameters  (s,  A,  7)  =  1, 0.95e'^°''*’^)} 

and  white  Gaussian  noise  in  an  8  x  8  snapshot  were  calculated  using  Monte  Carlo  simulations  each 
with  50  independent  experiments.  Several  2-D  methods  [10]  [25]  [38]  [59]  [2]  were  tested  for  accu¬ 
racy  without  interpolation,  with  focusing,  and  with  inverse  imaging.  These  methods  attain  either 
the  Cramer-Rao  bound  (CRB)  for  performing  the  estimation  independently  in  each  dimension  (a 
multi-trial  1-D  bound  or  the  1-D  by  1-D)  or  the  2-D  CRB  [11]  [52]  [66].  Figures  9  through  12  show 
the  Cramer-Rao  bounds  and  resultant  MSE  for  the  estimated  frequency  of  one  exponential.  Two 
solid  lines  are  displayed;  the  upper  is  the  2-D  Cramer-Rao  bound  and  the  lower  is  the  multi-trial 
1-D  bound. 

Degradation  away  from  the  Cramer-Rao  bound  is  observed  even  for  this  small  8x8  snapshot. 
The  ‘almost  rectangular’  assumption  fails  quickly.  When  the  sample  grid  is  greater  than  32  x 
32,  the  pixel  boundary  condition  is  violated  and  no  estimation  technique  performs  well  with  this 
assumption.  Each  of  the  interpolation  methods  improves  the  performance  of  the  algorithms  in 
spatial  frequency  estimation  (scatterer  location).  The  performance  of  some  estimation  methods, 
however,  is  degraded  and  no  longer  attains  the  Cramw-Rao  bound. 
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Mean  Squared  Error  of  arg(gamma)  mode  1 


Figure  9.  MSE  plot  for  almost  rectangular  assumption.  Techniques:  2-D  iterative  quadratic  max¬ 
imum  likelihood  (iqml2d),  matrix  enhancement  matrix  pencil  (memp),  2-D  method  of 
direction  estimation  (mode2d),  2-D  prony  (prony2d),  and  state  space  (statesp)  methods 


Mean  Squared  Error  of  arg(gamma)  mode  1 


Figure  10.  MSE  plot  for  scattering  center  focus.  Techniques:  2-D  iterative  quadratic  maximum 
likelihood  (iqml2d),  matrix  enhancement  matrix  pencil  (memp),  2-D  method  of  direc¬ 
tion  estimation  (mode2d),  2-D  prony  (prony2d),  and  state  space  (statesp)  methods 
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Mean  Squared  Error  of  arg(gamnna)  mode  1 


Figure  11.  MSE  plot  for  1-D  Lagrange  interpolation.  Techniques:  2-D  iterative  quadratic  max¬ 
imum  likelihood  (iqml2d),  matrix  enhancement  matrix  pencil  (memp),  2-D  method 
of  direction  estimation  (mode2d),  2-D  prony  (prony2d),  and  state  space  (statesp) 
methods 


Mean  Squared  Error  of  arg(gamma)  mode  1 


Figure  12.  MSE  plot  for  2-D  inverse  distance  interpolation.  Techniques:  2-D  iterative  quadratic 
maximum  likelihood  (iqml2d),  matrix  enhancement  matrix  pencil  (memp),  2-D 
method  of  direction  estimation  (mode2d),  2-D  prony  (prony2d),  and  state  space 
(statesp)  methods 
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The  benefit  of  combining  interpolation  and  estimation  is  easily  seen.  Not  only  will  the  inter¬ 
polation  method  then  be  appropriate  to  the  specific  method;  but  the  estimation  technique  could 
also  better  exploit  the  interpolation.  The  non-rectangular  data  grid  could  be  used  directly  by  the 
algorithm,  perhaps  to  form  a  projection  onto  the  signal  subspace.  Or  the  interpolation  results  such 
as  from  the  Lagrange  polynomial  could  be  used  instead  of  the  interpolated  data.  The  interpolation 
problem  is  not  as  significant  for  airborne  data  shown  in  Figure  3  where  the  k  —  w  domain  data  is 
evenly  sampled  in  the  direction,  however,  it  remains  to  be  seen  how  flight  path  deviations  and 
other  real-world  problems  effect  scatterer  location. 
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V.  Maximum  Likelihood  Estimation  of  1-D  Exponentials  in  Colored  Noise 

5.1  Overview 

Parametric  techniques  to  estimate  complex  exponential  signals  in  noise  play  an  important 
role  in  many  signal  processing  applications  such  as  direction  of  arrival  estimation  in  array  process¬ 
ing,  high  resolution  spectral  estimation,  and  radar  data  modeling.  An  overview  of  these  diverse 
techniques  is  found  in  [66] .  A  key  assumption  in  many  of  these  techniques  is  that  the  noise  spec¬ 
trum  is  flat  or  white.  When  the  noise  is  not  white,  more  accurate  estimation  of  the  exponentials 
is  possible.  White  noise  techniques  do  not  take  advantage  of  this  fact.  Additionally,  peaks  in  the 
noise  spectrum  may  be  incorrectly  identified  as  exponentials.  This  work  in  this  chapter  does  not 
rely  on  the  white  noise  assumption.  Instead,  exponential  estimation  techniques  are  based  on  the 
assumption  that  the  noise  is  unknown,  allowing  noise  with  spectral  variation  (colored  noise)  to 
be  parameterized  and  estimated.  Because  of  superior  statistical  performance,  maximum  likelihood 
(ML)  estimators  of  the  parameters  are  developed.  Several  approaches  [5]  [36]  [77],  [86]  to  find  the 
maximum  likelihood  estimate  of  exponentials  in  colored  noise  have  been  developed  that  model  the 
signal  as  random  amplitude  exponentials  (the  stochastic  model).  Most  of  these  methods,  however, 
require  several  instances  of  data  to  form  an  estimate  of  the  data  covariance  matrix.  Stochastic 
ML  techniques  are  examined  in  the  latter  half  of  this  chapter.  This  chapter  is  primarily  concerned 
with  single  data  instances.  In  this  case,  the  exponential  frequencies  and  amplitudes  can  be  con¬ 
sidered  unknown  constants  (the  deterministic  model).  An  approximation  to  the  ML  estimator  for 
the  parameters  of  a  colored  noise  deterministic  model  was  developed  in  [29].  Section  4.3  shows 
how  the  method  in  [29]  approximates  the  true  ML  estimator  derived  in  this  chapter.  This  chapter 
considers  the  estimation  of  the  most  general  case  of  exponential,  the  damped  exponential.  The 
methods  developed  can  be  easily  extended  to  undamped  exponentials  with  appropriate  constraints. 
The  methods  are  also  easily  extended  for  multiple  instances  of  data. 

This  chapter  employs  the  traditional  approach  to  finding  the  ML  estimators  of  the  exponential 
parameters  in  the  deterministic  model.  The  log-likelihood  function  of  the  data  and  parameters  is 
determined  by  taking  the  logarithm  of  the  probability  density  function  of  the  observed  data.  When 
a  parameter  can  be  decoupled  and  estimated  separately,  the  derivative  of  the  log-likelihood  function 
with  respect  to  the  parameter  is  determined.  The  derivative  is  then  set  equal  to  zero  in  order  to 
find  the  value  of  the  parameter  that  maximizes  the  likelihood  function,  i.e.,  the  ML  estimator  of 
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that  parameter.  ML  estimators  for  the  exponential  amplitudes  and  the  noise  variance  are  found 
in  this  manner.  When  these  estimates  are  substituted  back  into  the  log-likelihood  function,  the 
resulting  compressed  log-likelihood  function  is  a  nonlinear  function  of  the  exponential  frequencies 
and  the  noise  parameters.  Computationally  efficient  methods  to  find  the  parameters  that  maximize 
this  compressed  log-likelihood  function  are  developed  by  extending  the  white  noise  IQML  method 
[6]  [33]  [66].  Unknown  colored  noise  techniques  are  developed  for  the  cases  of  stationary  noise, 
autoregressive  (AR)  noise  (white  noise  driving  an  all-pole  filter),  and  autoregressive  moving-average 
(ARM A)  noise  (white  noise  driving  a  filter  with  poles  and  zeros). 

This  chapter  is  organized  as  follows.  Section  5.2  reviews  the  damped  exponential  model  and 
develops  the  maximum  likelihood  estimate  of  the  parameters  of  this  model,  develops  methods  for 
efficient  computation  of  the  deterministic  ML  estimates  for  several  colored  noise  models,  and  shows 
the  statistical  performance  of  these  new  methods  with  respect  to  the  Cramer-Rao  estimation  bound. 
Section  5.3  develops  efficient  methods  of  solving  the  stochastic  maximum  likelihood  colored  noise 
problem  and  shows  the  statistical  performance  of  these  methods. 

5.2  1-D  Deterministic  Maximum  Likelihood 

5.2.1  1-D  Damped  Exponential  Model.  This  section  describes  the  damped  exponential 

model  with  colored  noise.  Then  deterministic  maximum  likelihood  is  applied  to  estimate  the  pa¬ 
rameters  of  the  model.  The  estimation  of  the  damped  exponential  model  for  a  single  data  instance 
is  examined.  This  case  is  easily  modified  for  the  multiple  instance  case,  or  the  constrained  to  the 
unit  circle  case  by  techniques  described  in  [6]  and  [34]. 

The  complex  valued  data  samples,  t/[n]  for  n  =  0 . .  .iV  —  1,  of  the  damped  exponential  model 
are  given  by 

p 

2/W  =  n  =  0...N  -1,  (66) 

where  p  is  the  number  of  exponentials,  s,  is  the  amplitude  of  the  z**  exponential,  A”  = 
is  the  damped  exponential  with  damping  factor  a,-  and  exponential  frequency  /?,,  and  u;[n]  is  a 
zero-mean,  stationary  Gaussian  noise  sequence.  The  exponential  part  of  the  signal  is  considered  de¬ 
terministic  with  unknown  parameters,  {A  =[  Ai  A2  •••  Ap]^,s=[si  S2  •••  Sp  ]^}-  The 

randomness  of  the  data  is  due  entirely  to  the  normally  distributed  noise  sequence,  w((t^,?7)  ~ 


40 


N[0,  cr^-Rww(^)],  where  cr^  is  the  noise  variance,  and  i?ww(»?)  is  the  noise  covariance  matrix  with  q 
unknown  parameters  {77  =[771  772  •  •  ■  r)g  ]^}-  When  the  data  samples  are  arranged  in  a  column 

vector,  Equation  66  is  written  with  respect  to  the  model  parameters  9  —  {A,  s,cr^,  77}  as 

y{e)  =  G(\)s  +  w((T^  77),  (67) 

where  G(A)  is  an  iV  x  p  matrix  with  a  Vandermonde  column  of  each  of  the  damped  exponentials 
such  that 


-  [  gi(Ai)  g2(A2)  •  •  ■  gp(Ap)  ],  (68) 

gi(Ai)  =  [1  A.-  A?  •••  Af  f , 

and  s  is  a  column  vector  of  the  exponential  amplitudes.  Additionally,  w  is  complex  circularly 
symmetric  such  that  E{w*w}  =  R-wwiv)  and  E{w^w}  =  0,  where  (•)^  indicates  transpose 
and  (■)*  indicates  complex  conjugate  transpose.  Thus  the  data  y{9)  is  normally  distributed  as 
N[G(A)s,  (7277^^(77)]. 

5.2.2  Maximum  Likelihood  Estimators.  The  maximum  likelihood  estimate  of  the  param¬ 
eters  of  the  damped  exponential  and  colored  noise  model,  9  =  {A,  SjCr^,  77},  are  those  values  of  the 
parameters  that  maximize  the  probability  density  function  for  the  observed  data  sample  y.  For 
simplicity,  the  explicit  dependence  of  G  and  i?ww  on  the  model  parameters  A  and  77,  respectively, 
is  dropped  from  the  notation.  The  probability  density  to  be  maximized  is  the  circular  complex 
Gaussian  density, 

(69) 

|7r<7  ^ww| 

The  logarithm  is  a  monotonic  function,  thus  the  logarithm  of  the  probability  density  is  equivalently 
maximized  with  respect  to  parameters  in  A,  s,  <t^  and  77.  This  log-likelihood  function  for  the  density 
in  Equation  69  is 

C{G,  s,  0-2,  y)  =  -N  Inna^  -  In  |i7ww|  -  ^(y  -  Gs)*7?“^(y  -  Gs).  (70) 
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The  values  of  the  parameters  which  maximize  this  function  are  the  maximum  likelihood  estimates 
of  those  parameters.  The  derivatives  of  the  log-likelihood  function  with  respect  to  cr^  and  s  are 


^  +  ^(y  -  GsYR-Uy  -  Gs),  (71) 

^  =  -G*R-iGs  +  G*R-iy  (72) 

where  the  derivatives  with  respect  to  complex  vectors  are  defined  in  [74].  Setting  these  derivatives 
equal  to  zero  gives  the  maximum  likelihood  estimators  of  and  s, 

a^  =  j^(y-GsyR-i,{y-Gs),  (73) 

i  =iG*R:Y^G)-^G*R-ly.  (74) 

A  compressed  log-likelihood  function  can  then  be  constructed  by  substituting  the  ML  estimates  of 
and  s  into  the  log-likelihood  function, 

£,(^G,  RinYf'>y)  =  ^{G,  S,(T  TZwwjy)  1(7^=$^, s=s 

=  -N  In  ^y*(J  -  GiG*R-iG)-^G*RZiTR-Ul  “  G{G* R-^G^^G* R-i)y  -  In  |i?ww|  -  N. 

The  remaining  parameters  are  found  by  minimizing  the  negative  of  the  compressed  log-likelihood 
function.  Combining  terms  and  ignoring  constant  terms  the  minimization  is 

minln  |y*(7  -  G(G*fi:-i,G)-iG*72-i,)*i7-l,(/  -  G(G*i7-l,G)-iG*77-i,)yi?ww|.  (76) 

A,f) 

This  expression  is  quickly  simplified  when  we  define  a  projection  onto  the  signal  space 

Pl,g  =  L-*GiG*R-l,G)-^G*L-^  (77) 

where  the  positive  definite  covariance  matrix  is  decomposed  as  R-wvr  =  L*  L,  and  L~*  =  (L*)~^. 
The  minimization  is  then 

minln  |y*(7  -  L* Pl,gL-*)* R-i{I  -  L* PL,GL-*)yR^^\ 

X,ri 
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=  minln  |y*(r  i"*  -  L* Pl,gL-*)* R-UL* L-*  -  i*PL,GL"*)yi?wwl 

X,ft 

=  mm\n\y*  L-\l-PL,G)*LR-iL*  {I -PL,G)L-*yR„„\ 

X,ri 

=  minln  |y*i-i(/-PL.G)i'*yiZww|  (78) 

X,ri 

since  the  projection  {I—Pl,g)  is  idempotent  and  Hermitian  symmetric.  Next,  we  define  a  projection 
onto  the  space  orthogonal  to  the  signal  subspace  [9], 

Pl,a  =  LA(A*R„„A)-^A*L*,  (79) 

where  =  0.  The  iV  —  px  matrix  j4  is  formed  from  the  coefficients  (a  =[  aj  a2  •••  Op  V) 
of  a  polynomial, 

^(■^>)  =  ^  ^  +  •  ■  •  +  O'pX-  ^  =  0  (i  =  1 . .  .p),  (80) 

that  annihilates  all  of  the  exponential  columns  in  G 

dp  dp—\  '  *  •  Oq  0 

A*=  .  (81) 

0  dp  Up— 1  •  *  *  dQ 

The  Aj  i  =  1 . .  .p  are  simply  the  roots  of  the  polynomial  a(A).  Since  G  is  rank  p,  A  is  rank  N  —  p, 
and  G  and  A  are  orthogonal,  the  projections  Pt,G  and  Pl,a  form  a  resolution  of  the  identity  where 

Pl,g  +  Pl,a  =  Inxn-  (82) 

Thus,  the  minimization  in  Equation  78  in  terms  of  the  space  orthogonal  to  the  signal  space, 

minln  |y*i“^Pt,>ii“*yPww| 

a,r] 

=  minln  |y*A(j4*i?ww-d)“^A*yPww| 

a, I) 

=  ininln|y*yl(j4*i?„„j4)"^^*y/| +ln|Pwwl,  (83) 

will  give  the  ML  estimates  of  the  parameters  a,  the  annihilator  of  the  exponentials,  and  rj,  the  noise 
parameters.  This  is  a  nonlinear,  multidimensional  optimization  problem  which  could  be  solved  by 
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many  well  known  methods.  Notice,  however,  that  the  derivative  with  respect  to  a  of  the  terms  to  be 
minimized  in  Equation  83  does  not  include  the  second  term.  The  ML  estimator  of  the  exponentials 
does  not  depend  on  this  noise  related  term.  It  is  also  suggested  in  [30]  when  estimating  noise  that 
the  second  term  in  the  equation  simply  constrains  the  potential  solutions  (noise  poles)  away  from 
the  unit  circle. 

Performing  the  minimization  in  Equation  83  separately  with  respect  to  a,  the  minimization 

gives 

mmy*A{A*R„„A)-^A*y.  (84) 

a 

Further,  the  inverse  in  Equation  84  is  decomposed  as  (j4*72„vir^)~^  =  (A* L* LA)~^  —  (I/A)+(A*L*)+, 
where  (•)+  indicates  the  Moore-Penrose  pseudoinverse,  by  using  the  pseudoinverse  formula  (A*L*)+  = 
LA{A* L* LA)~^  and  (LA)'^  =  (A* L* LA)~^  A* L* .  The  minimization  is  then  stated  compactly  and 
efficiently  computed  as 


min  y*A(LA)+ (A*L*  )+A*  y 

a 

=  min||(^*L*)+^*y||f  (85) 

a 

2  ,  , 

where  {LA)* {LA)  is  the  Cholesky  decomposition  of  .4*/?wwj4  and  H-Hip  is  the  Frobenius  norm.  The 
following  section  shows  how  iiww  is  attained  and  how  the  minimization  in  Equation  85  can  be 
efficiently  computed. 

5.2.3  Conditional  Estimation  of  the  Exponentials.  This  section  introduces  a  new  method 
for  estimation  of  exponentials  in  unknown  colored  noise.  With  this  method  an  efficient  white  noise 
technique  is  extended  to  unknown  colored  noise  for  several  colored  noise  models.  Several  new 
algorithms  are  then  developed  to  perform  this  estimation. 

The  asymptotic  properties  of  maximum  likelihood  are  retained  when  Equation  84  is  opti¬ 
mized  with  respect  to  a  by  using  an  asymptotically  unbiased,  consistent  estimate  of  i2ww  Given  a 
consistent  estimate  of  i2ww,  R-w-w,  the  minimization  in  Equation  83  is  asymptotically  equivalent  to 

miny*yl(yl*.Rwwjd)“^A*y.  (86) 

a 
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As  is  the  case  with  other  ML  estimates,  in  addition  to  asymptotic  efficiency,  Equation  86  produces 
good  estimates  for  short  data  records  as  well.  For  the  case  of  white  noise,  IQML  [6]  [33]  [66]  is  a 
computationally  efficient,  iterative  method  to  find  the  parameters  a  in  Equation  86.  In  IQML,  the 
A*y  term  is  restated  as  Ya  where 


^*y 


0>p  •  ‘  ‘  CEq  0 


Oip  1  *  ‘  *  ^0 


Vp 

Vp-i  ■ 

yo 

yo 

yp+i 

Vp  ■ 

yi 

yi 

= 

yp 

VN-l 

VN-l 

yjv-p-i 

flo 

ai 


(87) 


With  Equation  87,  the  maximum  likelihood  minimization  in  Equation  86  becomes 


mina*Y*(A*.^wA)-iYa.  (88) 

a 

In  IQML,  A  is  constructed  from  the  elements  of  a  in  a  previous  iteration.  The  quadratic  form  in 
Equation  88  is  minimized  by  any  of  a  number  of  well-known  methods.  The  constraint  ||a||^  =  1 
guarantees  a  unique  solution  and  the  minimizing  a  is  then  the  eigenvector  associated  with  the 
minimum  eigenvalue  of  Y*(A*.Rwwj4)“^y. 

The  basic  iteration  used  in  all  the  algorithms  developed  in  this  section  is  summarized  at  the 
iteration  as 

mina;(Y*(A;(_iJ^wwAi:_i)“^Yai.  (89) 

The  kernel  y*(A’(_j.RwwAi:_i)“^Y  is  quickly  calculated  using  the  Cholesky  decomposition  of 
^fc-i-Rww^jfe-i  in  a  colored  noise  implementation  of  [23].  To  initialize  the  algorithm,  A  is  set 
equal  to  the  identity  matrix  I.  Thus,  the  first  iteration  whitens  the  data  with  the  covariance 
before  estimating  the  exponential  frequencies. 

To  demonstrate  the  achievable  performance  of  the  algorithms  developed  in  this  section  the 
known  colored  noise  version  of  IQML  where  Rww  =  Rww  in  Equation  89  was  implemented  and  is 
called  Colored  noise  IQML  (CIQML). 


Ya. 
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5.2.4  Estimation  in  Unknown  Colored  Noise.  Imposing  some  assumptions  in  order  to 
parameterize  the  noise,  allows  development  of  consistent  estimators  of  the  parameters.  In  all  the 
following  cases,  the  noise  is  assumed  ergodic.  Thus,  the  autocorrelation  sequence  of  the  noise  is 
consistently  estimated  [51]  as 

^  N—m 

^  ^  «^Mu;[n  +  m]  m  = -(A^  -  1) . .  .0  . . .  -  1.  (90) 

n=0 

The  problem  posed  by  Equation  89  is  how  to  initiate  the  iteration  with  a  good  estimate  of  a  and 
R-wvr.  Better  estimates  are  attained  when  the  signal  and  noise  are  estimated  simultaneously.  This, 
however,  is  not  always  possible. 

5.2. 4. i  Unknown  Colored  Noise  IQML  (UCIQML).  A  common  assumption  made 
is  that  the  noise  is  stationary.  This  is  the  case  solved  in  [77]  for  the  stochastic  model.  The 
matrix  iiww  is  then  Toeplitz,  Hermitian  symmetric  with  N  parameters.  These  N  parameters 
of  the  autocorrelation  function  are  based  on  N  real- valued  of  the  power  spectral  density  (PSD). 
Then,  together  with  the  undamped  exponential  signal,  4p-h  A  + 1  real- valued  parameters  are  to  be 
estimated.  For  the  single  data  instance,  this  problem  is  well-posed  since  N  complex- valued  samples 
of  data  are  available.  For  this  case,  the  signal  and  noise  are  estimated  separately  as  detailed  below. 

White  noise  exponential  estimators  are  consistent  in  colored  noise  [21].  That  is  asymptotically, 
the  residual  noise,  the  difference  between  the  estimate  and  the  observed  data,  can  not  be  used  to 
improve  the  exponential  estimate.  It  is,  however,  the  good  performance  of  ML  estimators  such  as 
IQML  for  short  data  records,  that  allows  us  to  estimate  the  signal  and  subtract  it  from  the  data 
to  estimate  the  noise  sequence.  From  Equation  67  the  noise  sequence  is  consistently  estimated  (see 
Appendix  B)  as 

w  =  y-Gs,  (91) 

where  G  contains  Vandermonde  columns  of  exponentials  estimated  by  IQML  and  s  is  estimated 
from  Equation  74  with  i?ww  =  /  as 


s  ={G*G)-^G*y  =  G+y, 


(92) 
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where  (•)■*■  indicates  the  pseudoinverse.  i2ww  is  simply  the  Toeplitz  covariance  matrix 

••• 

^ww[l]  i’ww[0] 

with  elements  consisting  of  the  autocorrelation  sequence  of  w,  for  =  —(N  —  1) ...  TV  —  1. 

With  the  estimate  of  a  from  IQML,  and  L  from  the  Cholesky  decomposition  of  i?ww,  the  iteration 
in  Equation  89  can  be  computed.  Additionally,  since  this  iteration  improves  the  estimate  of  a, 
and  consequently  the  estimate  of  s,  s  is  reestimated  at  each  iteration  by  Equation  74  and  w  is 
reestimated  by  Equation  91.  For  a  speedy  implementation,  the  autocorrelation  sequence,  fww[ru] 
for  m  =  —(TV  —  1) ..  .TV  —  1,  is  quickly  computed  with  the  Fast  Fourier  Transform  (FFT)  and 
estimates  are  quickly  computed  with  the  pseudoinverse  computed  by  the  QR  decomposition.  For 
example  Equation  74  becomes 

s=i{L*)+G)+iL*)+y.  (94) 

Just  as  in  IQML  of  the  iteration  in  Equation  89  is  not  guaranteed  to  converge,  however,  after 
a  small  number  of  iterations  (less  than  ten)  estimates  are  significantly  improved.  One  key  element 
of  sustaining  the  iteration  in  Equation  89  is  maintaining  the  positive  definiteness  of  the  kernel 
(A*Rww4)“^.  Thus,  only  positive  definite  estimates  of  the  noise  covariance,  J?ww)  are  used  in  the 
methods  developed  in  this  dissertation. 

5. 2. 4-2  The  UCIQML  Algorithm.  The  proceeding  algorithm  will  be  called  Unknown 
Colored  noise  IQML  (UCIQML)  and  is  summarized  as  follows: 

Step  1.  Estimate  the  signal  Gs  with  IQML.  Set  Aq  =  I  and  Rww  =  7  in  Equation  89.  Iterate 
several  times  to  attain  the  roots  of  a,  A,  for  i  =  1 ..  .p.  Form  the  matrix  G  from  the  A,-.  Then 
estimate  s  as 

s  =G'^y.  (95) 


fiww[0] 

fww[l] 

fiww[TV  -  1] 
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Step  2.  Estimate  the  noise  sequence  w  as 


w  =  y—Gs.  (96) 

Step  3.  Estimate  the  autocorrelation  sequence  of  the  noise  rww[^]  for  m  =  — (A^  — 1) . .  .N  —  l. 
Form  i^virw  and  its  Cholesky  decomposition  =  L*L. 

Step  4.  Estimate  new  A,-  for  i  =  I . .  .p.  Form  Ak-i  from  the  last  estimate  of  a,  and 
iterate  once 

min|l(A*_ir)+ya,l|;.  (97) 

Step  5.  Re-estimate  s  as 

s=((r)+G)+(r)+y.  (98) 

Step  6.  Repeat  Step  2  -  5  for  several  iterations,  or  until  estimates  do  not  significantly  change. 

5. 2. 4- 3  Other  Noise  Models.  In  addition  to  stationarity,  the  noise  may  be  further 
parameterized  as  the  stationary  output  of  a  filter  driven  by  a  white  noise  sequence.  Since  the 
available  data  is  limited,  an  appropriate  low-order  model  must  be  chosen  to  limit  the  number  of 
parameters  to  be  estimated.  The  principle  of  parsimony  [7]  suggests  that  to  better  estimate  the 
signal,  the  noise  be  modeled  with  as  few  parameters  as  possible. 

Two  well-known  models  that  provide  adjustable  low-order  noise  models  are  the  autoregressive 
(AR)  and  autoregressive  moving- average  (ARMA)  filter  models,  with  Power  Spectral  Densities 
(PSDs)  being  rational  functions  of  two  polynomials,  C{z)/ B{z).  In  the  next  sections,  specific 
solutions  for  AR  and  ARMA  noise  processes  are  explored.  It  is  emphasized  that  the  basic  colored 
noise  technique  exploited  in  this  research  is  very  general  and  any  appropriate  low-order  noise  model 
may  be  used  to  estimate  for  Equation  89. 

5. 2. 4- 4  Autoregressive  Unknown  Colored  Noise  IQML  (ARUCIQML).  The  follow¬ 
ing  solution  for  the  unknown  colored  noise  problem  assumes  an  AR  model  for  the  noise  (i.e.  the 
noise  contains  only  poles  and  no  zeros;  the  single  filter  coefficient  of  the  numerator  polynomial  is 
one).  If  the  noise  sequence  is  the  stationary  output  of  a  q*^  order  Infinite  Impulse  Response  (IIR) 
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filter,  then  the  noise  satisfies  an  autoregressive  (AR)  difference  equation 


hnw[n\  —  — 6ir(;[n  —  1]  ...  —  hqw\n  —  5]  +  e[u] 


(99) 


where  e[n]  is  a  white  noise  sequence.  The  AR  PSD  is 


p  (J<^\  _ 


(100) 


As  stated  earlier,  the  first  problem  posed  by  the  minimization  in  Equation  89  is  how  to 
attain  a  good  initial  estimate  of  Rww  c^nd  A  so  that  the  iteration  will  converge  to  the  true  A. 
It  is  well  known  that  either  AR  noise  poles  or  exponential  frequencies  may  be  estimated  with  the 
covariance  method  [29].  Additionally,  IQML  provides  consistent  ML  estimates  of  both  exponentials 
and  AR  processes  [66]  (see  Appendix  B).  The  only  difference  in  modeling  an  AR  noise  sequence  and 
exponentials  is  the  input  to  the  AR  filter  that  produces  them,  white  noise  or  an  impulse.  Figure 
13  shows  how  the  model  of  damped  exponentials  in  colored  noise  includes  both  these  cases.  If  we 
assume  that  the  colored  noise  can  be  described  with  an  autoregressive  (AR),  or  all-pole  filter  model 
with  known  model  order  q,  then  the  exponentials  and  noise  poles  may  be  estimated  simultaneously 
by  estimating  an  order  p  A  q  AR  polynomial  with  the  covariance  method  or  the  IQML  method. 


5. 2. 4-5  Approximating  the  Maximum  Likelihood  Optimization.  To  see  why  this  is 
true,  notice  that  Equation  84  may  be  generalized  as 


m^ny*AVFA*y,  (101) 

where  W  is  positive  definite.  A  generalization  to  any  positive  definite  matrix  W  is  possible  because 
the  solution  space  for  the  minimization  problem  is  defined  by  A*G  =  0  (the  null  space  A*).  The 
matrix  AW  A*  has  the  same  null  space.  Next,  decompose  the  inverse  covariance  matrix  as  in  [30] 
and  use  the  fact  that  is  centro-symmetric  (Hermitian  symmetric  and  symmetric  about  the 
anti-diagonal)  [45]  to  produce 

R-l,  =  =  iBP-^/^){p-^^^B*)  =  L-^L-*,  (102) 
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Figure  13.  Model  for  exponentials  in  ARMA  colored  noise 
where  B*  is  the  N  —  q  x  N  matrix 

1  6i  • • •  bg  0 

1  6l  •  •  •  bg 

1  •••  bill 

b\ 

0  1 

The  bi  are  the  AR  polynomial  coefficients,  b  =  [  1  bi  •  •  •  ] ,  the  for  i  =  1 . . .  fc  are  the 

order  AR  coefficients  from  the  Levinson  recursion  [30]  and 

P  =  diag([li  ...  ...  ^0  ]) 
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where  the  pk  foi  k  =  1 ..  .q  are  the  prediction  errors  calculated  at  the  iteration  of  the  Levinson 
algorithm.  The  following  approximations  of  Equation  84  asymptotically  produce  the  same  result 
and  approximate  the  maximum  likelihood  solution, 


miny*A(A*RvirwA)  ^A*y 

A 

mmy*A{A*{BP-^B*)-'^A)-'^A*y 

A 

(104) 

mmy*A(A*{KK*)+A)-^A*y  {q  <  p) 

A 

(105) 

miny*Ak{A*A)-^k*A*y 

A 

(106) 

mmy*Akk*A*y. 

(107) 

In  Equation  104,  except  for  the  last  q  rows  of  ,  both  B*  and  A*  perform  a  convolution 

on  the  data.  In  fact  since  pk+x  <  pk,  the  matrix  KK*  is  a  good  rank  N  —  q  approximation  of 
BP~^B* .  The  operations  performed  on  the  data  vector  y  are  now  clearly  seen  as  a  convolution 
(multiplication  by  j4*),  followed  by  a  deconvolution  (.A*)'*’,  followed  by  a  convolution  K* .  If  the 
order  of  these  operations  is  interchanged,  then  the  matrices  K*  and  (A*)"*"  may  be  replaced  by  the 
appropriately  sized  matrices  K*  and  (A*)"''  with  0{p/N)  error  (see  Appendix  A)  and  for  p  N, 
Equation  106  results.  This  is  the  ITEMS  minimization  [31]  (modified  for  damped  exponentials) 
with  iteration 

mmeilYKiAl_iAk-i)-^K*Y*eik.  (108) 

«)< 

Equation  107  indicates  that  the  covariance  method  will  also  provide  a  consistent  estimate  of  the 
exponentials  and  noise  poles.  Rewriting  the  minimization  in  terms  of  the  polynomials  whose  coef¬ 
ficients  are  contained  in  a  and  b  gives 


min||y(a*b)||^  (109) 

a*b 

where  *  indicates  the  convolution  of  the  polynomial  coefficients.  The  roots  of  the  polynomial  whose 
coefficients  are  a*b  are  the  exponentials.  A,-  for  i  =  1  ■  ■  p,  and  the  noise  poles,  Q  for  i  =  I ..  .q. 
The  difficulty  is  choosing  which  roots  are  exponentials  and  which  roots  are  noise.  The  amplitude 
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of  exponentials  is  estimated  from  Equation  74  with  -Rww  =  7  as 

s={H*H)-^H*y  =  H+y.  (110) 

If  H  contains  a  Vandermonde  column  for  each  root  of  the  order  p  +  q  AR  polynomial,  then  for 
high  signal-to-noise  ratios  (SNRs)  the  p  roots  with  the  largest  amplitude  will  he  the  exponentials. 
Equation  110  provides  a  good  indicator  of  signal  or  noise  down  to  near  zero  SNR. 

5. 2.4-6  Choosing  Signal  and  Noise  poles.  To  see  why  picking  roots  with  the  largest 
amplitude  is  a  reasonable  estimate  of  the  signal,  consider  the  maximum  likelihood  minimization 
for  white  noise 


mmly*  {l  -  G{G*G)-^G*)  y 

(Ill) 

max^y*G(G*G)-iG*y 

(112) 

max~y*GiG*G)-^G*G{G*G)-^G*y. 

(113) 

Substituting  back  in  the  ML  estimate  of  the  amplitude  gives  the  optimization 


max  —s*G*Gs. 
A  N 


(114) 


Thus  the  ML  estimate  maximizes  the  energy  of  the  estimated  signal.  Now,  consider  the  problem  of 
choosing  the  signal  and  noise  poles.  The  function  in  Equation  114  could  be  evaluated  at  all  (^p*) 
possible  signals  pole  combinations  and  the  combination  that  produced  the  maximum  value  chosen. 
A  similar  technique  could  also  be  employed  with  the  colored  noise  ML  optimization.  However, 
both  these  methods  are  computationally  intractable  for  large  p  or  q.  Thus,  a  suboptimal  approach 
is  needed.  Such  an  approximation  to  Equation  114  attained  by  ignoring  the  cross  terms  between 
signals  is 


max 
A  N 


1  P 


i|e 

gig*  Si, 


(115) 


t  =  l 


or  pick  the  p  poles  with  the  highest  energy.  For  undamped  exponentials  g*g  =  N,  and  Equation 
115  reduces  to 


p 


Si  Si, 


(116) 
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or  choose  the  poles  with  the  largest  amplitude.  For  damped  exponentials  the  appropriate  energy 
measure  is 

Note  that  while  these  approaches  are  suboptimal,  they  still  solve  an  optimization  problem  in  their 
own  right;  the  maximum  likelihood  estimate  of  a  single  exponential  signal.  Thus,  at  low  SNR  when 
the  poles  of  a  signal  consisting  of  closely  spaces  exponentials  have  merged  this  technique  provides 
better  estimates  of  the  merged  pole. 

Better  performance  at  low  SNR  can  also  be  attained  with  an  initial  estimate  of  the  p  q 
exponentials  and  noise  poles  from  IQML  as 

(118) 

Equation  109  also  suggests  that  a  comparable  performance  will  be  attained  by  increasing  the 
model  order  of  IQML  (IQML2)  to  model  both  signal  poles  and  noise  poles,  in  effect,  estimating  an 
AR(p  +  q)  process.  Results  show  that  this  is  true,  however,  picking  q  poles  as  noise  and  whitening 
the  data  with  Equation  89  still  achieves  better  performance.  That  is,  whitening  the  noise  poles 
over  all  frequencies  with  Equation  89  provides  better  performance  than  estimating  (eliminating 
from  the  signal)  the  subspace  of  the  noise’s  peak  frequencies  by  estimating  an  AR(p  4-  q)  process. 

5. 2. 4-7  The  ARUCIQML  Algorithm.  This  algorithm  will  be  called  Autoregressive 
Unknown  Colored  noise  IQML  (ARUCIQML)  and  is  summarized  as  follows: 

Step  1.  Estimate  p  exponentials  and  q  noise  poles  with  IQML.  Estimate  7,  for  i  =  1 . .  .p+  q 
and  form  the  N  x  p  +  q  matrix  H. 

Step  2.  Determine  which  are  exponentials  and  which  are  noise  poles.  Estimate  s  as 

s  =R+y.  (119) 

For  high  signal-to-noise  ratios  (SNRs)  the  p  exponentials  of  H  with  the  largest  amplitudes  are  the 
signal  exponentials,  Aj  for  i  =  1 . .  .p;  the  remaining  7,-  are  the  noise  poles,  Q  for  i  =  1 . .  .g. 

Step  3.  Construct  A  and  L.  Form  the  coefficients  of  a(A)  and  A  from  the  Aj-  for  f  =  1 . .  .p. 
Form  the  coefficients  of  b{Q  =  6o+^'iC+-  •  from  the  Q  for  i  =  1 . .  .q.  Form  the  autocorrelation 


mm 

a*b 


ik*A*)+Y(a-kh) 
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sequence  of  the  noise  from  the  AR  PSD.  Then  form  and  its  Cholesky  decomposition  = 
L*L. 

Step  4.  Estimate  new  A;  for  z  =  1 . .  .p.  Iterate  once 

rcdn\\(Al_,L*)+YB.k\\l.  (120) 

Step  5.  Form  A  from  the  last  estimate  of  a  and  repeat  Step  4  for  several  iterations,  or  until 
estimates  do  not  significantly  change. 

Step  6.  Estimate  s  as 

s  =  ((r)+G)+(r)+y.  (121) 

5.2. If. 8  Autoregressive  Moving  Average  Colored  Noise  IQML  (ARMAIQML).  The 

following  solution  for  the  unknown  colored  noise  problem  assumes  an  ARMA  model  for  the  noise 

(i.e.  the  noise  contains  poles  and  zeros).  If  the  noise  sequence  is  the  stationary  output  of  a  filter 
with  an  equal  number  of  poles  and  zeros  then  the  noise  satisfies  an  autoregressive  moving  average 
(ARMA)  difference  equation 

6ozn[n]  =  —biw[n  —  1]  ...  —  bqw[n  —  g]  (122) 

+Coe[n]  +  cie[n  -  1]  . . .  +  Cqe[n  -  q], 

where  e[n]  is  a  white  noise  sequence.  The  ARMA  PSD  is 

” - L.  (123) 

Damped  exponentials  are  a  deterministic  ARMA  process.  Thus,  IQML  is  again  used  to 
simultaneously  estimate  the  exponentials  and  the  noise  poles.  The  exponentials  and  noise  poles  are 
estimated  and  selected  by  amplitude  as  in  ARUCIQML.  ARUCIQML  also  estimates  the  exponential 
and  noise  pole  amplitudes  from  Equation  92.  However,  since  the  noise  amplitudes  are  stochastic,  the 
numerator  coefficients  of  the  ARMA  polynomial  can  not  be  determined  from  the  noise  amplitudes. 
Instead,  the  signal  and  noise  are  separated  as  in  UCIQML.  The  noise  is  filtered  (using  the  noise 
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poles  estimated  with  IQML)  to  attain  the  MA  part  (numerator)  of  the  noise  estimate  implied  by 
the  IQML  pole  estimate. 

As  in  ARUCIQML,  IQML  provides  consistent  estimates  of  the  exponentials  and  noise  poles 
and  approximates  the  maximum  likelihood  solution  (see  Appendix  B).  Similarly  to  [36]  the  ARM  A 
covariance  matrix  is  constructed  as 

Rarma  =  C*RarC  (124) 

where  Rar  isanA  +  pxA  +  p  AR  covariance  matrix  and  C  is  the  N  +px  N  convolution  matrix 
of  the  MA  (numerator)  polynomial, 


As  in  the  AR  case,  the  N  —  p  x  N  —  p  lower  triangular  Toeplitz  matrix  {K)'^  used  to  approximate 
R  —  {KK*)~^  commutes  with  C  and  A  (see  Appendix  A  and  B).  The  ARM  A  optimization  can 
thus  be  stated  as  miny* AKW K* A*y  with  a  positive  definite  weighting  matrix  W  and  IQML  will 
consistently  estimate  the  exponentials  and  noise  poles. 

First  for  the  ARMA  noise,  the  noise  and  signal  are  estimated  as  in  ARUCIQML.  Then,  the 
signal  and  noise  are  separated  as  in  UCIQML  as 

w  =  y-Gs  (126) 

where  G  contains  Vandermonde  columns  of  the  exponentials  estimated  by  ARUCIQML  and  s  is 
estimated  by 

s  =G+y.  (127) 
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The  noise  poles  were  also  estimated  and  the  AR  covariance  matrix  Rar  =  L*j^jiLar  is  constructed 
from  these.  The  estimated  noise  sequence  w  is  then  whitened  to  give  the  estimated  MA  noise 


sequence 

wma  =  iL*AR)'^y'-  (128) 

The  autocorrelation  sequence  of  the  M  A  coefficients  is  then  the  central  2q  +  l  values  of  the  autocor¬ 
relation  sequence  of  wma  ■  The  ARMA  covariance  matrix  is  constructed  from  autocorrelation 
sequence,  the  inverse  Fourier  transform  of  the  ARMA  PSD 

'p?  7.,,.  g-jwm 

P^u>ien  =  (129) 

With  the  estimate  of  a  from  the  signal  poles  and  L  from  the  Cholesky  decomposition  of  Rww, 
the  iteration  in  Equation  89  is  computed.  Again,  since  this  iteration  improves  the  estimate  of  a, 
and  consequently  the  estimate  of  s,  s  is  reestimated  at  each  iteration  by  Equation  127,  and  w  is 
reestimated  by  Equation  126  giving  a  new  MA  autocorrelation  sequence  vma  and  noise  covariance 
matrix  Rww 


5. 2.4-9  The  ARMAIQML  Algorithm.  This  algorithm  will  be  called  Autoregressive 
Moving  Average  colored  noise  IQML  (ARMAIQML)  and  is  summarized  as  follows: 

Step  1.  Estimate  p  exponentials  and  q  noise  poles  with  IQML.  Estimate  7j  for  i  =  1 . .  .p-h  g 
and  form  the  N  x  p  +  q  matrix  H. 

Step  2.  Determine  which  are  exponentials  and  which  7,-  are  noise  poles.  Estimate  s  as 

s  =H+y.  (130) 

For  high  signal-to-noise  ratios  (SNRs)  the  p  exponentials  of  H  with  the  largest  amplitudes  are  the 
signal  exponentials,  for  z  =  1 . .  .p.  The  remaining  ji  are  the  noise  poles,  Q  for  i  =  1 . .  .g. 

Step  3.  Construct  G,  A,  and  Ear-  Form  the  coefficients  of  a(A),  and  A  from  the  Aj  for 
i  —  I . .  .p.  Form  the  coefficients  of  b(Q  =  60  +  i*iC  +  •  ■  •  +  from  the  Q  for  i  =  1 . .  .q.  Form  the 
autocorrelation  sequence  of  the  AR  part  of  the  noise  from  the  AR  PSD.  Form  the  Toeplitz  matrix 
Rar  and  its  Cholesky  decomposition  Rar  =  L\hLar. 
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Step  4.  Estimate  the  MA  part  of  the  noise  sequence  wma-  Filter  out  the  AR  part  of  the  noise 


by 


WMA  =  iL*^R)+{y-Gs).  (131) 

Step  5.  Estimate  the  autocorrelation  sequence  of  the  MA  part  of  the  noise  fbr  = 

-(iV-  1)...A-  1. 

Step  6.  Estimate  the  autocorrelation  sequence  of  the  noise  rw„[m]  for 
m  =  —(N  —  1) ...  IV  —  1.  Retain  the  central  2g+ 1  values  of  the  cma  autocorrelation  sequence,  then 
calculate  the  autocorrelation  sequence  of  the  noise  from  the  Fourier  transform  of  the  ARMA 
PSD  in  Equation  129.  Form  and  its  Cholesky  decomposition  iZww  =  L*L. 

Step  7.  Estimate  new  A,-  for  i  =  I ..  .p.  Form  A  from  the  last  estimate  of  a  and  iterate  once 

™n||(A:_ir)+ya,||^.  (132) 

Step  8.  Estimate  s  as 

s  =((r)+G)+(r  )+y.  (133) 

Step  9.  Repeat  Step  4  -  8  for  several  iterations,  or  until  estimates  do  not  significantly  change. 

5.2.5  Results.  To  test  the  algorithms  in  previous  section,  instances  of  two  exponentials 
in  AR  and  ARMA  colored  noise  were  created.  For  the  first  example,  the  exponential  parameters 
were  (s.  A),  ={(1, 0.95e-^°  ®’^)i,  (1, 0.95e-’°  ®^’^)2}.  The  noise  was  a  two  pole  process  with  poles  ^  = 
{0.95e^°-^’^,  0.956-^ °  The  PSD  of  the  signal  and  noise  are  shown  in  Figure  14.  The  Mean  Square 

Error  (MSE)  in  estimating  the  exponentials  is  shown  in  Figure  15.  The  MSE’s  were  calculated  for 
length  N  =  S2  snapshots  using  Monte  Carlo  simulations  each  with  200  independent  experiments. 
All  of  the  algorithms  were  run  for  ten  iterations.  The  straight  solid  line  in  the  figures  is  the 
Cramer- Rao  estimation  Bound  (CRB)  which  is  described  in  Appendices  C  and  D. 

IQML,  IQML  with  double  the  known  model  order  (IQML2),  UCIQML  and  CIQML  were 
tested  with  this  data.  The  results  show  that  CIQML  follows  the  CRB.  The  MSE  of  UCIQML 
will  approached  the  CRB  asymptotically  as  the  length  of  the  sample  approaches  infinity.  For  the 
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Power  Spectrum  (0  dB  SNR) 


Figure  14.  PSD  of  two  exponentials  well  separated  in  frequency  from  AR  noise  poles  (Example  1) 


Mean  Squared  Error  of  arg(lambda)  mode  1 


Figure  15.  MSE  plot  for  one  of  two  exponentials  separated  from  AR  noise  (Deterministic  ML,  Ex¬ 
ample  1).  Techniques:  iterative  quadratic  maximum  likelihood  (iqml),  iqml  with  dou¬ 
ble  the  model  order  (iqml2),  known  colored  noise  iqml  (ciqml),  unknown  colored  noise 
iqml  (uciqml),  autoregressive  uciqml  (aruciqml),  and  iterative  estimation  of  mixed 
spectra  (items) 
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Figure  16.  PSD  of  two  exponentials  centered  near  in  AR  noise  poles  (Example  2) 

length  32  sample,  it  is  quite  close  to  the  CRB.  IQML  attains  the  CRB  for  white  noise  (i.e.  the 
CRB  for  colored  noise  is  higher  than  the  CRB  for  white  noise).  The  performance  of  both  IQML 
and  UCIQML  (because  it  is  based  on  the  signal  estimates  of  IQML)  breaks  down  at  a  moderate 
signal-to-noise  ratio  (30  dB).  The  performance  of  IQML2  does  not  experience  this  break-down  until 
much  lower  SNR.  This  suggests  the  techniques  exploited  in  ARUCIQML.  The  initial  estimate  for 
ARUCIQML  is  attained  in  the  same  manner  as  the  IQML2  estimate,  which  estimates  both  the 
signal  poles  and  the  noise  poles  with  its  increased  model  order.  The  ARUCIQML  algorithm  and 
the  approximate  ML  method  ITEMS  [29]  [31]  described  by  Equation  106  were  also  tested  on  this 
first  example.  The  exponentials  and  noise  poles  were  estimated  simultaneously  with  these  two 
algorithms  with  good  results.  Their  performance  matches  that  of  CIQML  and  tracks  the  CRB. 

For  the  second  AR  example,  the  exponential  parameters  were  changed  to  (s,  A)  ={(1,  0.95e-^°-^^”'), 
(1,0.956-'°  ^^’^)}.  These  exponentials  are  slightly  offset  in  frequency  from  the  noise  poles  at  ^  = 
(0.95e''°-^’^,  0.95e-'°  ‘^’^)  and  provide  a  more  difficult  estimation  problem.  The  PSD  of  this  example 
is  shown  in  Figure  16.  Results  are  shown  in  Figure  17.  In  this  example  the  increased  model  order 
of  IQML2  was  detrimental.  IQML2  performance  is  nonexistent  and  ITEMS  performance  is  poor. 
Using  ARUCIQML  in  both  AR  examples  results  in  improved  or  comparable  performance  with 
IQML. 
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Mean  Squared  Error  of  argflambda)  mode  1 


Figure  17.  MSE  plot  for  one  of  two  exponentials  centered  in  AR  noise  (Deterministic  ML,  Exam¬ 
ple  2).  Techniques:  iterative  quadratic  maximum  likelihood  (iqml),  iqml  with  double 
the  model  order  (iqml2),  known  colored  noise  iqml  (ciqml),  unknown  colored  noise 
iqml  (uciqml),  autoregressive  uciqml  (aruciqml),  and  iterative  estimation  of  mixed 
spectra  (items) 

For  the  ARMA  example,  the  exponential  parameters  were  (s.  A),-  ={(1, 0.95e'’°  ®’^)i,(l,  0.95e'^°  ®^’^)2}. 
The  noise  was  a  two-pole/two-zero  process  with  poles  ^  =  {0.95e'’°’^’^,  0.95e'^°  ‘*’'}  and  zeros  ^  = 
{0.95e'^°’^®’',  0.95e'^°  The  PSD  of  the  signal  and  noise  is  shown  in  Figure  18.  The  Mean  Square 
Error  (MSE)  in  estimating  the  exponentials  is  shown  in  Figure  19. 

Again,  CIQML  follows  the  CRB,  and  UCIQML  and  ARMAIQML  will  approach  the  CRB 
asymptotically.  IQML  attains  the  CRB  for  white  noise.  The  performance  of  both  IQML  and 
UCIQML  (because  it  is  based  on  the  signal  estimates  of  IQML)  breaks  down  at  a  moderate  signal- to- 
noise  ratio  (30  dB).  The  performance  (MSE“^)  of  the  unknown  colored  noise  algorithms  (UCIQML, 
ARUCIQML,  and  ARMAIQML)  is  bounded  above  by  the  performance  of  CIQML  and  below  by 
IQML  applied  to  the  colored  noise.  As  expected  when  the  model  used  better  matches  the  test  case 
(ARMAIQML  vs.  ARUCIQML,  ARUCIQML  vs.  UCIQML),  the  performance  is  better.  The  key 
benefit  of  UCIQML,  ARUCIQML,  and  ARMAIQML  is  that  the  better  (colored  noise)  CRB  that 
is  attained.  All  the  algorithms  derived  in  this  section  also  demonstrate  good  performance  in  white 
noise  with  exponentials  (s,  A)i  ={(1, 0.95e-’°  ®’^)i,(l,  0.95e''°  ®^’^)2}  as  shown  in  Figure  20. 
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Power  Spectrum  (0  dB  SNR) 


Figure  18.  PSD  of  two  exponentials  separated  from  ARMA  noise  (Example  3) 


Mean  Squared  Error  of  arg(lambda)  mode  1 


Figure  19.  MSE  plot  for  one  of  two  exponentials  in  ARMA  noise  (Examples  3).  Techniques:  itera¬ 
tive  quadratic  maximum  likelihood  (iqml),  known  colored  noise  iqml  (ciqml),  unknown 
colored  noise  iqml  (uciqml),  autoregressive  uciqml  (aru ciqml),  and  autoregressive  mov¬ 
ing  average  uciqml  (armaiqml) 
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Figure  20.  MSE  plot  one  of  two  exponentials  in  white  noise  (Example  4).  Techniques:  iterative 
quadratic  maximum  likelihood  (iqml),  unknown  colored  noise  iqml  (uciqml),  autore¬ 
gressive  uciqml  (aruciqml),  and  autoregressive  moving  average  uciqml  (armaiqml) 


5.B.6  Conclusions.  Computationally  efficient  maximum  likelihood  estimators  for  colored 
noise  have  been  developed.  These  estimators  retain  all  the  asymptotic  properties  of  maximum 
likelihood  and  produce  good  results  for  short  data  records  as  well.  The  ARUCIQML  algorithm 
attains  the  CRB  and  is  efficient  for  short  data  records.  The  UCIQML,  an  order  MA  colored 
noise  estimator,  and  ARMAIQML  algorithms  do  not  attain  the  CRB.  This  is  expected  since  the 
colored  noise  estimates  are  based  on  finite  order  sample  covariance  estimates  which  are  not  efficient 
for  MA  and  ARM  A  processes  [71]. 


5.3  Stochastic  Maximum  Likelihood 

This  section  introduces  a  new  technique  to  estimate  exponentials  in  colored  noise  based  on 
stochastic  maximum  likelihood.  This  section  simplifies  the  solution  for  the  colored  noise  stochastic 
model  given  in  [86]  from  a  global  search  to  an  iterative  quadratic  solution.  This  simplification  uses 
the  projection  onto  the  orthogonal  subspace  developed  in  the  previous  section  and  the  methodology 
of  the  white  noise  stochastic  ML  method.  Method  of  Direction  Estimation  (MODE)  [69] . 

The  results  of  [86]  allow  known  and  unknown  colored  noise  IQML-like  algorithms  for  the 
maximum  likelihood  solution  of  the  stochastic  model  to  be  derived.  The  data  in  the  stochastic 


62 


model  is  distributed  as 


y  ~  JV(0,  Ryy) 

Ryy  =  GRssG*  +  (T^Ry,y,  (134) 

The  general  result  of  [86]  finds  the  maximum  likelihood  estimate  of  the  unknown  parameters 
with  the  minimization  with  respect  to  the  parameters  G  and  -Rww  of 

£(G,  y)  =  Min  \PRyyP*  +  0-^(1  -  P)Ry,y,\  +  MN  (135) 

where  P  =  L* Pl^gL~*  and  =  (A^  —  p)“^tr(/  —  P)RyyR~\,.  Observe  the  term  inside  the  deter¬ 
minate  is  the  estimated  covariance  matrix  separated  into  signal  subspace  and  noise  subspace  terms. 
In  terms  of  the  deterministic  signal,  y  =  Gs  -|-  <tL*€  and  Ryy  =  yy*  the  minimization  is 

min  In  |(Gs  -b  <TG{G*R-iG)-^G*L-h){Gs  +  (TG(G*R-iG)-^G*  L-^)*  +  &\l-  P)R^„\.  (136) 

LtfG 

When  the  signal  and  noise  are  uncorrelated  this  gives 

minln  |Gss*G*  +  (t^G{G* R-^G^'^G* L'^e* L-*G{G* R-iG)-^G*  +  or\l  -  P)Ry,y,\.  (137) 

Then  taking  the  expected  value,  E{  },  of  Equation  137  gives 

minln  |(Gss*G*  -b  a^PR^^  +  -  F)i7ww|.  (138) 

L,G 

The  first  term  in  this  minimization,  the  signal,  is  fixed,  thus,  the  minimization  is  satisfied  when 
is  minimized.  Observe  that  is  still  minimized  when  the  orthogonal  projection  Pl,a  is  used  in 
place  of  Pg,a-  Replacing  Pg,a  with  Pl,a  iu  Equation  135  and  some  simplification  gives 

minln  \L*PL,AL-*RyyiL*PL,AL-*y  +  d-^{I  -  rPL.AT“*)Rww| 

L,A 

=  minln  |Rww^(j4*Rww^)~^^*Ryy-4(-4*Rww-^)~^-d*Rww 
L,A 

+a\l  -  R^^A{A*R^^A)-^A*)R^^\ 

=  minln  |J?ww^(jd*Rww^)~^^*Ryy^(.d*i?ww-'4)“^A*Rww 
L,A 

-&^R„„A{A*R^^A)~^A*R^y^A{A*R„y^A)~'^A*R^^  -b  d^Rwwl 
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=  minln  |i?ww^(^*-Rww^)  ^A*{Ryy-&^R.„„)A{A*R„„A)  +  ^^-Rwwl 

L,A 

=  minln  |i?ww-A(A*i2ww^)~^^*(-Ryy  —  a-^ R-wv/)A{A* RTf,-„A)~^ A*  +  a-'^I\  +  In  |i?ww| 

L,A 

=  minln \A{A* R^^A)-'^A*{Ryy  -  d-^R„„)  +  +  In  {\AB  +  I\  =  \BA  +  I\) 

LyA 

=  minln \{A*L*)+A*{Ryy  -  R„„)A{LA)+  +  a'^I]  +  In  |i?ww|  (139) 

L,A 

Notice  again  that  the  second  term  in  Equation  139  does  not  depend  on  A.  Thus  this  term  can  be 
eliminated  from  the  minimization  to  find  A.  Since  the  logarithm  is  a  monotonic  function  the  mini¬ 
mization  in  Equation  139  is  equivalent  to  one  with  the  logarithm  removed.  Then,  for  a  consistent 
estimate  of  iiww,  Rww,  the  minimization  is 

min \{A*L*)+A*{Ryy  -  a^R„„)A{LA)+  +  a^I\  (140) 

L,A 

Note  the  quantity  in  the  determinant  is  positive  semidefinite.  Then,  the  inequality 

\R\  <  (^r  (141) 

for  R  positive  semidefinite  [41]  indicates  that  the  minimization  in  Equation  140  is  equivalent  to 

mmtr  (j^A*L*)'^A*{Ryy  —  Ryry^)A[LA)'^  +  (142) 

For  white  noise  this  is  equivalent  to  the  weighted  subspace  minimization  given  in  [50] .  It  is  shown 
in  [71]  that  the  stochastic  ML  model  may  be  applied  with  consistency  to  deterministic  data.  In  the 
form  of  Equation  142,  the  effect  of  the  minimization  on  the  signal,  y  —  Gs  +  (TL*e,  is  to  minimize 

E{tr  ((A*L*)+A*(Ayy  -  &'^R^^)A{LAy 

=  E{tr  {{A*L*)+A*Gss*G*A{LA)+  +  a^{A* L*)+ A* L*€€* LA{LA)+  -  &^I  +  &^l)} 

=  tv{(T'^{A*L*)+A*L*E{ee*}LA{LA)+) 

=  o-2tr(Pi,^) 

=  {N-p)a‘^  (143) 
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where  Eww  =  yy*,  tr(j4B)  =  tr(5j4),  and  A*G  —  Q.  Thus,  the  variance  of  the  noise  is  minimized. 
Observe  that  the  difference  from  the  data  model  is  the  way  in  which  the  effect  of  the  noise  on 
the  exponential  signal  estimation  (choice  of  minimizing  A)  is  minimized.  In  the  deterministic 
model,  a  whitening  filter  minimizes  the  noise.  In  the  stochastic  model,  an  estimate  of  the  noise  is 
subtracted  from  the  signal  prior  to  estimation  of  the  minimizing  A.  This  is  readily  accomplished  on 
the  whitened  data  in  an  extension  of  the  white  noise  technique  MODE  [69]  to  known  colored  noise 
which  will  be  called  Colored  noise  MODE  (CMODE).  MODE  is  described  in  Appendix  E.  Including 
our  previously  canceled  L*L~*  and  L~^L,  and  taking  the  principal  components  (pc)  of  the  whitened 
data  (the  eigenvectors  associated  with  the  p  largest  eigenvalues)  gives  the  minimization 

inintr  {(A* L*)+A* L*pc\^L-*{Ryy  -  LA{LA)+'j  (144) 

or 

rnintr  (^A{A*R„„A)-'^A*L*pc  {l~*(fiyy  -  (145) 

where  cr^  may  be  quickly  estimated  as  the  average  of  the  remaining  TV  —p  eigenvalues.  Notice,  that 
since  is  now  considered  any  consistent  estimate  of  cr^,  and  because  of  the  additive  property  of  the 
trace,  the  final  term  in  Equation  142  may  be  dropped.  This,  in  turn,  implies  that  any  weighting 
of  may  be  subtracted  in  Equation  145  so  long  as  the  matrix  R  in  Equation  141  remains  positive 
semidefinite.  A  discussion  of  the  optimum  weighting  is  included  in  the  next  section. 

Equation  145  is  minimized  in  a  similar  fashion  to  IQML  by  creating  a  data  matrix  Y  in 
Equation  87  to  form  the  minimization  of  a  quadratic  in  a.  Here,  also,  is  the  first  compromise 
necessary  for  a  single  data  instance.  The  full  covariance  matrix  Ryy  is  not  available;  yy*  is  rank  1 
and  a  principal  components  analysis  is  meaningless.  A  smaller  {N—pxN—p)  rank-p+1  covariance 
matrix  Ryy  =  YY*  may  be  estimated  from  overlapping  segments  of  the  data  with  Y  from  Equation 
87.  In  1-D  this  only  increases  the  number  of  samples  required.  In  2-D  estimation  this  will  force  a 
compromise  between  true  2-D  exponential  estimation  and  1-D  by  1-D  estimation. 

5.3.1  Estimation  in  Unknown,  AR  and  ARMA  Colored  noise.  Stochastic  algorithms 
corresponding  to  each  of  the  deterministic  algorithms  were  developed  and  tested  on  the  same  set 
of  examples  as  in  Section  4.4.  The  stochastic  algorithms  are  built  around  the  iteration  in  Equation 
145.  An  eigendecomposition  of  the  whitened  outer  product  of  the  data  matrix  in  Equation  87, 
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(£*)■'■  RL'^  =  {L*)^YY* L'^  —  VDV* ,  gives  p  signal  eigenvectors  L*VpW^I‘^  where  Vp  contains  the 
p  columns  of  V  associated  with  the  p  largest  eigenvalues  in  D,  and  W  is  the  optimal  weighting  of  the 
eigenvectors  [71]  (see  Appendix  E).  The  p  signal  eigenvectors  are  used  in  place  of  the  data  matrix  of 
the  colored  noise  version  of  IQML.  The  best  performing  and  most  computationally  efficient  method 
of  finding  principal  components  involved  taking  the  SVD  of  {L*)'^Y  .  This  quickly  identifies  the 
eigenvectors  Vp  and  the  square  roots  of  the  eigenvalues  in  D.  Other  aspects  of  the  algorithms  (how 
the  noise  covariance  is  consistently  estimated)  are  the  same  as  the  deterministic  versions. 

The  optimum  weight  matrix  W  for  a  large  number  of  data  instances  is  [50] 

W  =  {D-  =  {D-  or‘^I){I  -  (146) 

The  effect  of  this  weighting  matrix  can  be  understood  in  terms  of  the  inequality  in  Equation  141 
which  is  derived  from  the  dominance  of  the  arithmetic  mean  of  the  positive  and  zero  eigenvalues 
of  the  matrix  R  in  Equation  141  over  the  geometric  mean  of  these  eigenvalues.  This  inequality 
becomes  an  equality  when  the  eigenvalues  are  equal.  Then,  as  — »•  D  (low  SNR)  the  weighting 

term  in  Equation  146  goes  to  zero  and  the  approximation  of  using  the  trace  in  Equation  142  becomes 
more  accurate. 

The  same  two  mode  test  cases  as  in  Section  4.4  were  run  with  MODE,  and  its  known  colored 
noise  (CMODE)  and  unknown  colored  noise  variants  (UCMODE,  ARUCMODE,  ARMAMODE). 
The  results  are  shown  in  Figures  21  through  24.  When  multiple  instances  of  the  data  are  avail¬ 
able  and  when  the  exponentials  in  the  signals  are  correlated  these  methods  can  show  improved 
performance  over  the  IQML  based  methods  [71].  For  the  single  data  instance,  however,  these  al¬ 
gorithms  do  not  in  general  outperform  the  IQML  beised  methods.  In  some  instances,  the  MODE 
based  algorithms  do  show  better  performance  (UCMODE).  As  in  most  instances  of  the  determin¬ 
istic  algorithms,  the  unknown  noise  versions  of  the  algorithms  very  closely  approach  the  known 
colored  noise  version’s  performance.  The  MODE  based  algorithms,  however,  all  fall  slightly  below 
the  CRB.  This  is  due  the  fact  that  an  extra  p  samples  are  used  to  effectively  estimate  the  signal 
subspace,  thus  the  optimization  in  Equation  145  is  based  on  a  set  of  length  N  —  p  data  samples 
while  the  deterministic  optimization  is  performed  on  length  N  data  sample. 
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Mean  Squared  Error  of  arg(lambda)  mode  1 


Figure  21.  MSE  plot  for  two  exponentials  widely  separated  from  AR  noise  (Stochastic  ML,  Ex¬ 
ample  1).  Techniques:  method  of  direction  estimation  (mode),  known  colored  noise 
mode  (cmode),  unknown  colored  noise  mode  (ucmode),  and  autoregressive  ucmode 
(arucmode) 


Mean  Squared  Error  of  arg(lambda)  mode  f 


Figure  22.  MSE  plot  for  two  exponentials  centered  in  AR  noise  (Stochastic  ML,  Example 
2).  Techniques:  method  of  direction  estimation  (mode),  known  colored  noise 
mode  (cmode),  unknown  colored  noise  mode  (ucmode),  and  autoregressive  ucmode 
(arucmode) 
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Mean  Squared  Error  of  arg(lambda)  mode  1 


Figure  23.  MSE  plot  for  two  exponentials  in  ARMA  noise  (Stochastic  ML,  Example  3).  Tech¬ 
niques;  method  of  direction  estimation  (mode),  known  colored  noise  mode  (cmode), 
unknown  colored  noise  mode  (ucmode),  autoregressive  ucmode  (arucmode),  and  au¬ 
toregressive  moving  average  ucmode  (armamode) 


Mean  Squared  Error  of  arg(lambda)  mode  1 


Figure  24.  MSE  plot  for  two  exponentials  in  white  noise  (Stochastic  ML,  Example  4).  Tech¬ 
niques:  method  of  direction  estimation  (mode),  unknown  colored  noise  mode  (uc¬ 
mode),  autoregressive  ucmode  (arucmode),  and  autoregressive  moving  average  ucmode 
(armamode) 
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VI.  Estimation  of  2-D  Exponentials  in  Colored  Noise 

6.1  Overview 

This  chapter  develops  extensions  to  two  dimensions  of  the  one-dimensional  maximum  likeli¬ 
hood  exponential  estimation  techniques  developed  in  the  previous  chapter.  This  extension  is  not 
completely  straight  forward  because  of  the  unique  nature  of  two-dimensional  (2-D)  signals.  Two- 
dimensional  discrete  time  signal  processing  theory  has  been  well  developed  [14] ,  and  many  of  the 
characteristics  follow  directly  from  one-dimensional  discrete  time  signal  processing.  However,  one 
dimensional  spectral  estimation  techniques  that  involve  polynomial  rooting  fail  to  be  directly  ap¬ 
plicable  due  to  the  absence  of  a  fundamental  theorem  of  algebra  in  two  dimensions.  The  roots  of 
a  polynomial  in  two  dimensions  are  not  unique  and  efficient  root  finding  techniques  do  not  exist. 
This  fact  has  lead  to  the  development  of  several  novel  methods  of  2-D  spectral  estimation.  This 
section  briefly  reviews  the  2-D  techniques  that  have  been  developed  in  order  to  understand  how 
the  two-dimensional  maximum  likelihood  solution  effectively  exploits  both  dimensions.  First,  to 
discuss  the  2-D  techniques,  the  following  signal  model  of  a  single  instance  of  data  with  exponentials 
in  two  dimensions  is  used.  Two-dimensional  damped  exponential  data,  y,  in  noise,  n,  containing  p 
exponentials  may  be  described  as 

p 

y[mi,m2]  =  +  n[mi,m2]  0  <  mi  <  Mi  -  1  0  <  m2  <  M2  -  1  (147) 

Z  =  1 

where  the  Aj-  and  71  are  the  complex  exponential  frequencies  and  the  s,-  are  the  complex  amplitudes 
of  the  exponentials.  Arranging  the  data  y  in  an  Mi  x  M2  matrix  Y  gives 

Y  =  GSH'^  Y  N,  (148) 

where 

^=[gi  g2  •••  gp  ]  gi  =  ^Mi(Ai),  (149) 

S  =  diag([  Si  S2  •  •  •  Sp  ])  (150) 

H  =  [hi  h2  •••  hp  ]  hi  =  'JM,(Ti),  (151) 

and 

^m(z)  =  [1  ^  ^2  ...  (152) 
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AT  is  a  random  field  of  noise,  and  S  is  diagonal.  Another  less  restrictive  model  of  the  data  is  also 
possible  if  all  the  elements  of  the  matrix  S  are  populated.  In  this  case,  the  model  is  separable 
and  1-D  techniques  can  be  efficiently  applied  in  each  dimension.  In  this  case,  exponentials  with 
any  combination  of  spatial  frequencies  in  each  dimension  are  possible  (a  1-D  by  1-D  model).  This 
model  then  requires  a  strict  pairing  of  spatial  frequencies  if  two-dimensional  exponentials  at  unique 
frequencies  are  present  in  the  data.  Conversely,  the  restricted  (full  2-D)  model  where  S  is  diagonal 
does  not  require  pairing  since  its  solution  is  unique. 

6.2  Two-Dimensional  Techniques 

Several  two-dimensional  exponential  estimation  techniques  are  introduced  here  that  follow  a 
1-D  by  1-D,  or  a  2-D  model.  A  detailed  analysis  of  the  performance  of  these  techniques  is  contained 
in  [52].  The  techniques  discussed  here  are  those  that  do  not  involve  search  of  the  spatial  frequency 
plane,  but  instead  find  the  spectral  peaks  through  some  more  efficient  means. 

The  first  technique  from  [2]  is  developed  from  a  state  space  representation  of  the  data  and 
uses  a  singular  value  decomposition  of  the  matrix 

y  =  c/sy* 

where  U  and  V  are  unitary  matrixes  spanning  the  column  space  and  row  space  of  Y  respectively, 
and  E  is  a  diagonal  matrix  of  the  singular  values  of  Y.  The  matrices  G,  S,  and  H  are  estimated 
from  U,  E,and  V,  respectively.  The  spatial  frequencies  are  estimated  separately  in  each  dimension 
from  U  and  V  by  exploiting  the  rotational  invariance  of  exponentials  [57] .  As  in  all  the  remaining 
techniques  requiring  pairing,  the  frequencies  are  efficiently  paired  [52]  by  row  and  column  of  the  p 
largest  amplitudes  in  the  least  squares  estimate  of  the  matrix 

5  =  G+y(F^)+.  (153) 

This  method  attains  the  multi-trial  1-D  CRB  given  in  [10]. 

The  second  technique.  Matrix  Enhanced  Matrix  Pencil  (MEMP)  from  [25]  also  exploits  ro¬ 
tational  invariance  to  find  the  spatial  frequencies.  This  method  also  requires  pairing,  however,  it 
efficiently  estimates  the  2-D  exponentials  of  the  restricted  2-D  model  by  simultaneously  using  the 
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data  from  all  the  columns  or  rows,  and  attains  the  2-D  CRB  [10].  In  MEMP,  the  frequencies  are 
estimated  from  a  block  Hankel  matrix  consisting  of  blocks  which  are  the  Hankel  data  matrix  of  the 
1-D  covariance  method  [30]  for  each  column  of  data  .  This  block  matrix  is  in  effect  an  estimate  of 
the  Cholesky  decomposition  of  the  full  Mi  M2  x  Mi  M2  covariance  matrix  of  the  data,  which  has 
been  ordered  columnwise  into  one  vector.  However,  since  only  one  data  instance  is  available  and 
some  samples  are  used  up  to  make  the  matrix  rank  p,  the  block  matrix  is  not  the  full  size  2-D 
covariance  matrix. 

The  2D  Prony  method  [59]  is  an  extension  of  the  1-D  Prony  technique  [30].  The  covariance 
method  is  applied  simultaneously  to  all  the  columns  (rows)  of  data,  in  effect,  considering  each 
column  (row)  as  a  new  data  instance.  This  method  uses  the  one-dimensional  technique  of  over¬ 
modeling  (estimating  more  frequencies  than  the  data  contains)  and  selects  the  two-dimensional 
frequency  combinations  that  produce  the  highest  energy  signal. 

The  2D  MODE  technique  [38]  estimates  the  full  2-D  covariance  matrix  from  the  outer  product 
of  a  column-ordered  vector  of  the  data  (a  rank  one  estimate  for  the  single  data  instance).  Separately 
averaging  of  the  elements  of  the  covariance  matrix  associated  with  one  column  (row)  of  S  allows  full 
rank  covariance  matrices  for  each  dimension  to  be  used  in  the  1-D  stochastic  maximum  likelihood 
technique,  MODE  [69].  In  a  later  section,  the  2D  MODE  technique  is  extended  to  create  colored 
noise  techniques  based  on  stochastic  maximum  likelihood. 

For  high  signal-to-noise  ratios  (SNRs),  all  of  these  techniques  perform  well  with  respect  to 
the  appropriate  CRB.  Of  these  diverse  techniques,  2D  MODE  [38]  and  2D  IQML  [10]  (discussed  in 
detail  later)  can  be  related  to  a  maximum  likelihood  solution  of  the  two-dimensional  exponential 
estimation  problem  in  colored  noise,  and  only  MEMP  [25]  and  2D  IQML  [10]  achieve  the  2-D  CRB 
[52]. 

6.3  Deterministic  Maximum  Likelihood 

In  this  section,  the  two-dimensional  maximum  likelihood  problem  is  formulated,  then,  via  this 
formulation  the  techniques  are  extended  to  the  unknown  colored  noise  problem.  The  2-D  maximum 
likelihood  solution  follows  almost  directly  from  1-D  solution  when  the  columns  of  the  data  matrix 
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are  stacked  (this  is  the  tiectorize  operation).  Equation  148  then  becomes 


y  =vec(y)  =  (H  ®  G)s  +  n 


(154) 


or 

y=Fs  +  n  y  ~  N(f  s,cr2i?„„)  (155) 

where  ®  denotes  the  Kronecker  product,  s  =  vec(5),  and  n  =  vec(lV).  Equation  148  implies  that 
the  2-D  harmonics  lie  on  a  1-D  by  1-D  grid.  Removing  this  restriction  by  forcing  S  to  be  diagonal 
eliminates  the  elements  of  s  in  equation  154  where  Sj  =  0.  The  new  Mi  M2  x  p  matrix  F  becomes 

•f’  =  [hi®gi  h2(g)g2  hp(g)gp]  (156) 

where  s  now  contains  just  the  diagonal  elements  of  S.  Additionally,  the  noise  is  circular  Gaussian 
such  that  E{n*n}  =  cr^Rnn  and  E{n^n}  =  0,  where  (•)*  denotes  complex  conjugate  transpose. 

Analogous  to  the  1-D  case,  the  probability  density  and  log-likelihood  functions  are 


(157) 


and 


£(E,s,o-^Rnn;y)  =  -MiM2ln7r<7-2  -ln|R„n| - ^(y  -  Es)*R„^(y  -  Fs).  (158) 


The  zeros  of  and  are  the  maximum  likelihood  (ML)  estimates  of  and  s. 


d2  = 


1 


Ml  M2 


(y-FsTR-^iy-Fs) 


s  ={F*R-^F)-^F*R-^y  =  L*  PL,FL-*y 


(159) 

(160) 


where  the  positive  definite  covariance  matrix  is  decomposed  as  i2„„  =  L*  L  and 
Pl,f  =  L~* F{F* R~^F)~^ F* where  L~*  —  {L*)~^.  Substituting  Equations  159  and  160  into 
Equation  158  gives  the  compressed  log  likelihood  function 


jO{F,  Rnn;  y)  =  -Ml  M2  In  -  L*Pl,fL-*)*R-^{I  -  L*PL,FL-*)y  -  In  -  Mi  M2. 

JVL\M.2 


(161) 
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The  maximum  likelihood  estimates  of  the  remaining  parameters  are  then  found  by  maximizing 
C{G,  Rnn',  y)  or  equivalently  minimizing 


minln|y*L  ^(I  -  Pl,f)L  *yRnn\ 

=  minln|y*Z,“^PL,»vT“*y-Rnn|  (162) 

where  the  orthogonal  projection  to  Pl,f  is 

Pl,w  =  LW{W*RnnW)-^W*L*  (163) 

where  the  column  space  of  W  spans  the  subspace  orthogonal  to  the  signal  subspace,  W*F  —  0 
and  W,  or  its  eigenspace  equivalent  is  full  column  rank.  Then  using  the  Cholesky  decomposition  of 
Rnn  =  L*L,  and  full  column  rank  and  full  row  rank  definitions  of  the  Moore-Penrose  pseudoinverse, 


A+  =  {A*A)-'^A*  and  yl+  =  A*{AA*)-^  respectively, 

Pl,w  =  LW{W*RnnW)-'^W*RnnWiW*RnnW)-^W*L*  (164) 

=  LW{W*RnnW)-^W*L*LW{W*RnnW)-^W*L*  (165) 

=  LW{LW)+{W*L*yW*L*  (166) 

=  LW{W*Rn^W)+W*L*.  (167) 


Thus,  the  projection  can  be  constructed  from  W,  even  if  W  is  not  full  row  or  column  rank,  as  is 
the  case  in  the  following  section. 

Again  as  in  the  1-D  case  the  minimization  in  Equation  162  may  be  written  as 

minln  \y* Pl,w L~*y\  +  In  |i?„n|  (168) 

where  parameters  in  W  can  be  maximized  independently  of  the  ln|i?nn|  term.  Thus,  given  a 
consistent  estimate  of  ,  Rn^  the  parameters  in  W  that  minimize  Equation  168  are  found  from 

imny*WiW*^nW)+W*y.  (169) 
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6.3.1  Known  Colored  Noise  Maximum  Likelihood.  A  matrix  W  that  spans  the  subspace 
orthogonal  to  the  signal  subspace  is  derived  in  [9]  as  the  ((2Mi  —  1)(M2  —  p)  +  Mi  —  l)x  Mi  M2 
matrix 


Im,  ®  c; 

[Imi  [0]  ®  •D*_i  —  [0|/mi]  <8>  i 


(170) 


where  7  =  [  0i  O2  •••  Op_i  ]  and  the  (M2-p)xM2  matrixC*,  and  (M2-P+I)  XM2 

matrix  77*  _i  are  defined  from 


Cl 


Ck  ■■■  Cl  Co  0 


0  Ck  ■■■Cl  Co 


(171) 


and 


dk-i 


di  do  0 


D 


*  _ 

ib-1  — 


(172) 


[  0  dk-i  ■  ■  ■  di  do  J 

The  upper  block  of  W*  annihilates  the  harmonics  in  one  direction  A,- ,  and  the  lower  block  of  W* 
predicts  the  harmonics  in  the  second  dimension  jj ,  based  on  the  harmonics  annihilated  in  the  first 
dimension.  The  minimization  is  performed  with  IQML-like  iterations, 


mmf*Y*{w:_iRnnWi-i)+Yfi,  (173) 

•t 

where  f  =[  ,  c  =  [  cq  ci  •••  CkV,^=[do  di  •••  dji,_i  and  Y  is  defined 

such  that 

Y{=W*y.  (174) 

The  pseudoinverse  in  Equation  173  requires  a  significant  amount  of  memory.  An  new  equivalent 
implementation  based  on  the  relation 


{w*wy  =  {iww*)+wy{ww*)+w 


(175) 
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developed  in  Appendix  F,  requires  only  one  quarter  the  memory.  Then,  the  minimization  in 


Equation  173  becomes 


min 

ti 


i{LWi-iW*_  1  L*yLWi- 1  )yfi 


2 

F 


(176) 


Computations  can  be  reduced  further  from  0{M^)  (M  =  Mi  =  M2;  p  M)  for  the  SVD  in 
(W*W)'^  when  this  is  replaced  with  the  normal  inverse  by  using  a  permutation  matrix.  The 
matrix  P  chooses  Mi  M2  —  p  linearly  independent  rows  of  the  banded  matrix  W*  [12].  With  the 
(M1M2-P)  X  (Ml M2  —p)  matrix  PW*WP*  Equation  173  becomes 


mmf;+,Y*P*iPW:RnnWiPT"PH+i.  (177) 

fi+i 

Each  of  the  Kronecker  products  in  the  definition  of  W*  is  a  upper  band  matrix  with  the  last  having 
the  largest  bandwidth  M.  Thus  there  exists  a  matrix  PW*  with  bandwidth  M  and  a  band  Cholesky 
solution  to  Equation  177  is  O(M^)  [19]. 

6.S.2  2-D  Unknown  Colored  Noise  Maximum  Likelihood.  This  section  develops  a  new 
technique  to  estimate  2-D  exponentials  in  unknown  colored  noise  by  extending  the  2D  IQML 
technique.  When  the  noise  is  white,  =  I,  the  algorithm  to  perform  the  minimization  in 
Equation  173  is  called  2D  IQML.  This  algorithm  achieves  the  2D  CRB  [9]  [11]  [52].  As  in  the  1-D 
case,  due  to  the  Wold  decomposition  the  2-D  CRB  for  the  exponential  parameters  is  the  same  for 
unknown  colored  noise  and  known  colored  noise.  Also  as  in  the  1-D  case,  /?„„  can  be  replaced  with 
a  consistent  estimate.  Ran,  and  the  asymptotic  properties  of  maximum  likelihood  minimization  in 
Equation  168  will  be  retained.  The  unknown  colored  noise  variants  for  2-D  now  follow  directly 
form  the  1-D  cases. 

If  the  noise  sequence  is  assumed  stationary  then  the  noise  sequence  is 


n  =  y—Fs 


(178) 


where  F  is  estimated  by  2D  IQML  and  s  is  estimated  by 

s  ={F*F)-^F*y  =  F+y.  (179) 
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Then  for  stationary  noise,  i?„„  is  the  block  Toeplitz  covariance  matrix  [28]  whose  elements  are  the 
autocorrelation  coefficients  of  N  in  Equation  148.  N  is  attained  by  reversing  the  vec  operation  and 
extracting  length  Mi  segments  of  the  vector  n  and  ordering  them  into  a  matrix.  The  autocorrelation 
coefficients  of  N  are  then  the  2Mi  —  lx  2M2  —  1  matrix  Vnn-  The  matrix  f^N  is  quickly  attained 
by  taking  the  magnitude  squared  of  the  two-dimensional  FFT  of  N  zero-padded  to  2Mi  -  1  by 
2M2  —  1.  Then  the  2-D  autocorrelation  matrix  is  the  block  Toeplitz  matrix 


— 

where 


i?o  R—  1 

R-M2 

Ri  Rq 

Rm^  Rm^—i 

Ro 

»^ArAr[-l,i] 

rNN[-Mi,i] 

fwjv[0,i] 

rNN[-{Mi  -  l),i] 

(180) 


(181) 


[  fiv7v[Mi,f]  fivAf[Mi  -  l,i]  fjvjv[0,i]  J 

With  the  estimate  of  f  from  2D  IQML  and  L  from  i?„„  the  iteration  in  Equation  173  can  be 
computed.  Additionally,  since  this  iteration  improves  the  estimate  of  s,  and  consequently  the 
estimate  n  at  each  iteration  s  is  estimated  by  Equation  160  and  n  is  reestimated  by  Equation  178. 


6.3.2. 1  The  2D  UCIQML  algorithm.  This  techniques  will  be  called  2D  UCIQML 
and  is  summarized  as  follows: 

Step  1.  Estimate  the  signal  Fs.  Attain  the  two  dimensional  frequencies  (A,  7),-  for  i  =  1 . .  .p 
from  2D  IQML.  Form  the  matrix  F  from  the  frequencies.  Then  estimate  s  as 


s  =E+y. 


(182) 


Step  2.  Estimate  the  noise  sequence  n  as 


n  =  y-Es. 


(183) 
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Step  3.  Esiimaie  the  autocorrelation  of  the  noise  rjvAT-  Form  the  block  Toeplitz  matrix  Rnn 


and  its  Cholesky  decomposition  Ran  =  L* L. 


Step  4.  Estimate  new  (A, 7),-  for  i  =  1 . .  .p.  Form  Wi-i  from  the  last  estimate  f,  fi_i,  and 
iterate  once 


min 

fi 


{{LWi-iWUL*)+ LWi-i)Yfi 


2 

F 


(184) 


Step  5.  Re-estimate  s  as 

s=((r)+F’)+(r)+y.  (185) 


Step  6.  Repeat  Step  2  -  5  for  several  iterations,  or  until  estimates  do  not  significantly  change. 


6.3.3  2-D  Exponential  Noise.  This  section  develops  a  new  spectral  model  for  2-D  signals. 
Algorithms  to  exploit  an  AR  or  ARMA  model  in  two  dimensions  do  not  follow  as  readily  as  in  one 
dimension  since  these  models  do  not  have  the  same  properties  in  two  dimensions  as  they  do  in  one 
dimension.  The  rational  polynomial  model  for  colored  noise  is  less  appealing  due  to  the  absence 
of  a  fundamental  theorem  of  algebra  in  two  dimensions.  The  2-D  AR  model  is  not  efficient  on  a 
finite  region  of  support  and  the  set  of  AR  coefficients  describing  a  specific  PSD  is  not  unique.  The 
estimated  parameters  of  the  rational  polynomial  noise  models  in  two  dimensions  vary  according  to 
the  technique  used  to  estimate  the  model  parameters. 

Instead,  a  suitable  low-order  noise  model  with  unique  parameters  is  a  filter,  whose  impulse 
response  is  the  sum  of  damped  exponentials.  In  1-D,  this  model  is  a  subset  of  the  ARMA  noise 
model.  The  sum  of  q  damped  exponentials  is  synonymous  with  the  1-D  ARMA(g  —  l,q)  model. 
The  one-dimensional  ARMA(g  —  1,  g)  model  transfer  function  (z-transform  domain)  in  polynomial, 
pole-zero  and  partial  fraction  form  may  be  written  as 

c{z)  nL~/(i-o^-^)  ^  s. 

Biz)  -  eLo  hz-'^  -  ro=i(i  -  ^  (1  -  ^  ^ 


The  1-D  AR(g)  case,  excepting  the  constant  term,  the  coefficients  of  numerator  polynomial  in  1-D 
ARMA(g  —1,5)  model  equal  zero  and 


Cjz)  _  _ 1 _ _ _ 1 _ _  Sh 

Biz)  -  EU  hz-’^  ~  nLi(l  -  PkZ-^)  "  S  (1  -  ■ 


(187) 
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In  both  cases,  the  impulse  response  of  the  transfer  function  H{z)  =  C(^)/5(z)  is  sum  of  q  expo¬ 
nentials, 

9 

/i[n]  =  ^  sjbjSj  n  =  0...Ar  — 1.  (188) 

k-l 

Thus,  the  filter  for  an  AR  process  may  be  implemented  with  either  the  IIR  filter  coefficients,  bk 
iov  k  =  I ..  .q,  OT  the  FIR  filter  coefficients,  c/  =  Sk^l  ior  I  —  0  ... N  —  1.  The  exponential 

amplitudes,  sj,  for  =  1 . .  .9,  can  be  adjusted  to  create  any  ARMA(5  —  1,  q)  filter. 

If  this  idea  is  expanded  to  two  dimensions,  then,  the  pole  locations  of  any  2-D  filter  can  be 
fixed  and  the  2-D  AR  coefficients  are  not  needed.  Thus  the  impulse  response 

9 

h[mi,  m2]  =  0  <  mi  <  Ml  -  1  0  <  m2  <  M2  -  1  (189) 

k  =  l 

represents  the  coefficients  of  a  2-D  filter  with  known  pole  locations.  This  filter  is  the  unique  pairing 
of  two  1-D  ARMA  spectra  and  suitably  constrains  our  noise  model.  Equation  189  represents  a  new 
2-D  spectral  model.  For  the  case  where  g  =  1  this  defines  a  2-D  AR  filter.  In  1-D,  a  higher  order 
AR  filter  could  be  formed  by  cascading  (convolving)  several  AR(1)  filters.  In  2-D,  however,  this 
results  in  a  filter  with  poles  at  all  combinations  of  the  spatial  frequencies  in  each  dimension.  That 
is,  the  filter  becomes  separable  and  less  general.  In  the  separable  noise  model,  the  2-D  covariance 
matrix  is  independent  in  each  dimension  and  is  the  Kronecker  product  of  the  covariance  matrix  for 
each  dimension. 

Run  =  Rm,  <S>  RMi  ■  (190) 

Colored  noise  techniques  developed  in  this  chapter  will  be  tested  with  both  the  exponential  and 
separable  noise  model. 

6.3.4  Exponeniia!  Noise  Maximum  Likelihood.  This  section  develops  another  new 

technique  to  estimate  2-D  exponentials  in  unknown  colored  noise  by  extending  the  2D  IQML 
technique.  With  the  Exponential  noise  model,  the  2-D  extension  of  the  ARUCIQML  technique 
now  follows  directly.  Again  the  only  difference  in  modeling  the  exponential  noise  sequence  and 
exponentials  is  the  input  to  the  filter  that  produces  them,  white  noise  or  an  impulse.  Thus  the 
exponentials  and  noise  poles  may  be  estimated  simultaneously  by  estimating  p  +  q  exponentials 
with  the  2D  IQML  method.  To  choose  which  exponentials  are  signal  and  which  are  noise,  the 
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exponential  amplitudes  are  estimated  from  Equation  179.  The  p  exponentials  with  the  largest 
amplitude  are  chosen  to  represent  the  signal  and  the  remaining  exponentials  are  chosen  as  the 
noise  poles.  The  2-D  autocorrelation  matrix,  J?nn)  for  exponential  noise  can  then  be  formed.  With 
the  estimate  of  f  from  the  p  2-D  exponentials  chosen  as  signal  with  2D  IQML,  and  L  from  the 
Cholesky  decomposition  of  Ran,  the  iteration  in  Equation  173  can  be  computed. 


6. 3.4.1  The  2D  ARIQML  algorithm.  This  algorithm  will  be  called  2D  ARIQML 
because  of  its  similarity  to  the  1-D  method  ARUCIQML  and  is  summarized  as  follows: 

Step  1.  Estimate  p  exponentials  and  q  noise  poles  and  their  amplitudes.  Attain  the  two 
dimensional  frequencies  and  amplitudes  (A,  7,  s),-  for  i  =  1 . .  .p  -f  g  from  2D  IQML 

Step  2.  Determine  which  (X,j,s)i  for  i  =  l...p  +  q  are  exponentials  and  which  are  noise 
poles.  For  high  signal-to-noise  ratios  (SNRs)  the  p  exponentials  of  with  the  largest  amplitude  are 
the  signal  exponentials  (A,  7,  s),-  for  i  =  1 . .  .p;  the  remaining  (A,  7,  for  f  =  p  -|-  1 . . . g  are  the 
noise  poles. 

Step  3.  Construct  W  and  L.  Form  W  from  the  (A, 7),-  for  i  =  1 . .  .p.  Calculate  the  autocorre¬ 
lation  coefficients  of  exponential  noise  poles.  Then  form  the  Toeplitz  matrix  Ran  and  its  Cholesky 
decomposition  R„„  =  L*L. 

Step  4.  Estimate  new  (A,  7)^  for  i  =  1 . .  .p.  Iterate  once 


min 

fi 


(191) 


Step  5.  Form  W,_i  from  the  last  estimate  of  f,  fi_i,and  repeat  Step  4  for  several  iterations, 
or  until  estimates  do  not  significantly  change. 

Step  6.  Estimate  s  as 

s  =  ((L*)+E)+(L*)+y.  (192) 


6.4  Stochastic  Maximum  Likelihood 

This  section  extends  the  2D  MODE  technique  to  colored  noise  and  develops  two  new  tech¬ 
niques  to  estimate  2-D  exponentials  in  unknown  colored  noise.  2-D  stochastic  maximum  likelihood 
techniques  hinge  on  an  accurate  estimate  of  the  2-D  data  covariance  matrix.  For  a  single  data 
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instance,  the  outer  product  estimate  of  the  covariance  matrix  is  only  rank  one  and  thus  of  limited 
utility.  Other  estimates  of  the  2-D  covariance  matrix  such  as  the  outer  product  of  the  MEMP  data 
matrix  are  not  full  size  and  thus  provide  limited  accuracy.  The  assumption  of  stationarity  or  of 
separability  leads  to  a  full  rank  and  full  size  2-D  covariance  matrix  (Equation  180  and  190,  respec¬ 
tively).  But,  because  the  2-D  covariance  matrix  has  2M1M2  —  Mi  —  M2  + 1  independent  parameters 
[14],  and  only  Mi  M2  samples  are  available  in  a  single  data  instance,  the  full  2-D  stationary  and 
separable  covariance  matrix  is  of  limited  use  for  high  resolution  estimation  of  exponentials.  The 
independent  components  of  the  covariance  matrix  in  Equation  190,  however,  provide  the  basis  for 
computationally  less  intensive  colored  noise  algorithms  based  on  a  1-D  by  1-D  formulation. 

The  2-D  techniques  based  on  2D  IQML  in  the  previous  section  are  computationally  intensive, 
©(MfMl),  especially  as  Mi  and  M2  increase.  A  computational  savings  is  attained  when  colored 
noise  techniques  are  based  the  separable  model  such  as  the  stochastic  maximum  likelihood  technique 
2D  MODE  [38].  The  2D  MODE  technique  minimizes  an  alternative  form  of  Equation  169 

imntr  (^Pl^wL~* RyyL~^^  (193) 

where  Pl,w  is  given  in  Equation  167  and  for  the  single  data  instance  Ryy  =  yy*.  The  2D  MODE 
technique  assumes  that  the  exponentials  are  any  possible  pairing  of  the  estimated  frequencies  (a 
separable  model).  Thus  the  matrix  F  is  defined  as  the  Kronecker  product  of  G  and  H, 

F  =  H®G.  (194) 

To  extend  the  2D  MODE  technique  to  colored  noise  we  will  also  assume  that  the  noise  is  separable 
as  in  Equation  190  then 

Pl,w  =  I  —  Pl,f 

=  I-L-*F{F*R-^F)-^F*L-^ 

=  I  -  {L-^\  ®  L-^\){H  ®  G){{H*  ®  G*){Rm\  ®  R-J;){H  ®  G)r\H*  ®  G*){Ll,\  ®  L-Jj 

=  I  -  Lm*H{H*RiIh)-^H*LI}^  ®  L-^\G{G*R-^\G)-^G*L-^\ 

=  I-Ph®Pg  (195) 
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where  Pg  =  and  Pq  =  LM*GiG*RM\G)-^G*LM\.  Then  expressing 

Ph  and  Pq  in  terms  of  their  orthogonal  projections  Pp  and  Pq  . 

Pl,w  =  I  -Ph®Pg 

=  I-{I-P^)®{I-P^) 

=  P^®  J  +  /®Pi--P^®Pi'  (196) 

where  I  is  an  identity  matrix  appropriately  sized  for  each  Kronecker  product.  Equation  193  is  then 

mintr((P^®7  +  7®P^-P^®Pi-)i“*PL-^)  (197) 

where  R  =  yy*.  2D  MODE  considers  the  Kronecker  product  P^  ®  Pg  to  be  higher  order  and 
it  is  neglected.  Thus,  the  minimization  to  find  the  and  the  7,-  can  be  performed  separately. 
Effectively,  the  trace  of  the  blocks  of  the  2-D  covariance  matrix  are  averaged  across  one  dimension. 
Since  both  the  2-D  covariance  matrix  estimated  by  the  outer  product  of  the  vectorized  data,  and 
the  2-D  covariance  matrix  of  a  stationary  noise  process  have  special  structure  (an  outer  product, 
and  Toeplitz  blocks,  respectively),  the  average  across  blocks  in  each  dimension  can  be  computed 
directly  without  the  computationally  intensive  and  memory  consuming  task  of  forming  the  2-D 
covariance  matrices.  The  minimization  for  Aj  is 

inintr  {{I ®  P^)L-* RR-^)  =  inintr  {P^Lm\RxLm\)  (198) 

where  Rx  =  Pmimi  and  Ru  =  R{kM2  +  l :  {k  +l)M2,  IM2  +  I  '■  (I  +1)M2)  is  an  M2  x  M2 

submatrix  of  R.  The  minimization  for  j,  is 

mintr  ((P^  ®  I)L~*RL~^)  =  mintr  [P^L]^  RyL]^  )  (199) 

Ij  7j 

where  (TTy)*;  =  tr(Pji;;).  Both  minimizations  are  then  accomplished  by  the  1-D  colored  noise 
MODE  techniques  discussed  in  the  previous  chapter.  The  1-D  techniques  CMODE,  UCMODE, 
and  ARUCMODE  are  used  almost  directly  in  the  2-D  techniques  2D  CMODE,  2D  UCMODE,  and 
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2D  ARMODE.  These  minimizations  then  have  the  form 


mina*tr*(A*_iR„„Ai_i)-i[/ai,  (200) 

at 

where  Rm  =  L*j^Lm  and  U  =  Llf[  ■  ■  ■  Up  V  i®  defined  such  that 

Ula  =  A*Ui.  (201) 

Here  [/  =  [  U2  •  •  •  Up  ]  contains  the  p  principal  eigenvectors  of  or 

The  components  of  the  noise  covariance  matrix  L,  Lmi  and  Lms,  are  estimated  by  forming  R„„ 
from  Equation  180  as  in  2D  UCIQML  and  2D  ARIQML,  then  averaging  in  the  same  manner 
that  LJ^^RxL]^^  and  are  formed  from  L~* RL~^ .  The  harmonics,  Aj  or  7j,  are  then 

the  roots  of  the  polynomial  given  by  the  coefficients  in  a.  Pairing  the  Ai  and  7j  is  done  with 
Equation  153.  The  SVD  of  L'^^RxL'^^  and  dominates  the  computations  making  this 

an  0{M^)  technique  where  M  =  Mi=M2. 

6. 4-0-2  The  2D  UCMODE  algorithm.  The  stochastic  maximum  likelihood  technique 
based  on  unknown  colored  noise  will  be  called  2D  UCMODE  and  is  summarized  as  follows: 

Step  1.  Estimate  the  signal  Fs.  Attain  the  two  dimensional  frequencies  (A,  7),  for  z  =  1 . .  .p 
from  2D  MODE.  Pair  the  frequencies  by  least  square  amplitude  as  in  Equation  153,  however,  allow 
no  frequency  to  be  used  more  than  once  (to  attain  distinct  frequencies).  Form  the  matrix  F  from 
the  frequencies.  Then  estimate  s  as 

s  =E+y.  (202) 

Step  2.  Estimate  the  noise  sequence  n  as 

n  =  y-Fs.  (203) 

Step  3.  Estimate  the  autocorrelation  of  the  noise  rpfi<f.  Form  the  Toeplitz  matrices  Rmi  ^ind 
Rm2  )  S'Wd  find  their  Cholesky  decompositions,  Lm^  and  Lms  - 
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Step  4.  Estimate  new  (A,  7),-  for  i  =  1 . .  .p.  Form  vlj_i  from  the  last  estimate  a,  ai_i,  for 
each  dimension 

min  ^  (204) 

Bi  F 

Step  5.  Re-estimate  s  as 

s=iiL*)+F)+{L*)+y.  (205) 

Step  6.  Repeat  Step  2-5  for  several  iterations,  or  until  estimates  do  not  significantly  change. 

6. 4-0. 3  The  2D  ARMODE  algorithm.  The  stochastic  maximum  likelihood  tech¬ 
niques  based  on  exponential  colored  noise  will  be  called  2D  ARMODE  and  is  summarized  as 
follows: 

Step  1.  Estimate  p  exponentials  and  q  noise  poles  and  their  amplitudes.  Attain  the  two 
dimensional  frequencies  and  amplitudes  (A,  7,5),-  for  z  —  1 . .  .p  +  q  from  2D  MODE.  Pair  the 
frequencies  by  least  square  amplitude,  however,  allow  no  frequency  to  be  used  more  than  once 
(distinct  frequencies). 

Step  2.  Determine  which  (A,7,  s)j  for  z  =  \ . .  .p  +  q  are  exponentials  and  which  are  noise 
poles.  For  high  signal-to-noise  ratios  (SNRs)  the  p  exponentials  of  with  the  largest  amplitude  are 
the  signal  exponentials  (A,  7,  s),-  for  z  =  1 . .  .p;  the  remaining  (A,  7,  s)<  for  z  =  p  -|-  1 . . .  g  are  the 
noise  poles. 

Step  3.  Construct  A  and  L.  Form  A  from  the  (A,  7),-  for  z  =  1 . .  .p.  Form  the  autocorrelation 
coeflBcients  of  the  exponential  noise  poles.  Form  the  Toeplitz  matrices  Rm,  and  Rms,  and  find  their 
Cholesky  decompositions,  Lmj  and  Lms- 

Step  4.  Estimate  new  (A,  7),-  for  z'  =  1 . .  .p.  Form  Aj_i  from  the  last  estimate  a,  ai-i,  for 
each  dimension 

min  {A*_iL*M)+Uai  ^  (206) 

a;  F 

Step  5.  Form  Ai_i  from  the  last  estimate  of  a,  aj_iand  repeat  Step  4  for  several  iterations, 
or  until  estimates  do  not  significantly  change. 

Step  6.  Estimate  s  as 

s  =  ((T*)+F)+(r)+y.  (207) 
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Power  Spectral  Density 


Figure  25.  PSD  of  two  exponentials  in  exponential  noise 


6.5  Results 

Simulations  were  run  with  all  techniques  in  this  chapter  on  both  the  exponential  and  separable 
noise  models.  In  the  simulations,  the  2-D  colored  noise  techniques  estimated  two  modes  (s,  A,  7),'  = 

(e'^°'^’^,0.95e'’°  ‘*’^,  1)2}  in  exponential  noise  and  separable  expo¬ 
nential  noise  with  two  frequencies  in  each  dimension  {s„,Xn,jn)i  =  {(1, 0.97e“'’°  0.97e"'^‘''^’^)i, 

(1, 0.97e-^°’^’^,  0.97e-^°  '^’^)2}.  The  PSD  of  the  signal  in  exponential  noise  is  shown  in  Figure  25  and 
the  signal  in  separable  noise  is  shown  in  Figure  26. 

Data  instances  of  8  x  8  samples  were  used.  The  results  for  the  deterministic  techniques  are 
shown  in  Figure  27  and  28.  The  results  for  the  stochastic  techniques  are  shown  in  Figure  29  and 
30.  Mean  Square  Error’s  (MSE’s)  were  calculated  using  Monte  Carlo  simulations  each  with  200 
independent  experiments.  All  of  the  algorithms  were  run  for  three  iterations.  The  straight  solid 
line  in  the  figures  is  the  deterministic  2D  CRB  which  is  described  in  [12]. 

The  performance  (MSE“^)  of  2D  CIQML  shows  that  the  2-D  deterministic  techniques  can 
attain  the  CRB.  However,  no  efficient  estimator  of  the  2-D  noise  was  found  and  the  colored  noise 
techniques  attain  the  CRB  only  asymptotically.  The  unknown  colored  noise  technique  2D  UCIQML, 
which  is  not  efficient  in  1-D  (see  Chapter  4),  and  the  pole  estimating  technique  2D  ARIQML,  which 
may  not  be  efficient  for  a  finite  number  of  samples  in  2-D,  did  not  attain  the  2-D  CRB  [30].  The 
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Power  Spectr<y  Density 


Figure  26.  PSD  of  two  exponentials  in  separable  noise 


Mean  Squared  Error  of  arg(gamma)  mode  2 


Figure  27.  Estimation  error  for  one  of  two  exponentials  in  exponential  noise  (Deterministic  ML). 

Techniques:  2-D  iterative  quadratic  maximum  likelihood  (iqml2d),  known  colored 
noise  iqml2d  (ciqml2d),  unknown  colored  noise  iqml2d  (uciqml2d),  and  exponential 
noise  uciqml2d  (ariqml2d) 
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Mean  Squared  Error  of  arg(gamma)  mode  2 


Figure  28.  Estimation  error  for  one  of  two  exponentials  in  separable  noise  (Deterministic  ML). 

Techniques:  2-D  iterative  quadratic  maximum  likelihood,  fast  version  (iqml2f),  known 
colored  noise  iqml2f  (ciqml2f),  unknown  colored  noise  iqml2f  (uciqml2f),  and  expo¬ 
nential  noise  uciqml2f  (ariqml2f) 

techniques  did  improve  performance  over  the  white  noise  techniques  2D  IQML  and  2D  MODE. 
An  exception  is  the  2D  ARIQML  technique  which  does  not  estimate  the  noise  poles  accurately 
enough  in  comparison  to  the  signal  poles.  Thus,  the  inaccurate  noise  covariance  estimate  degrades 
2D  ARIQML  estimation  performance.  This  is  true  despite  the  fact  that  in  estimates  of  the  noise 
only  2D  ARIQML  estimated  the  noise  better  than  2D  ARMODE,  which  shows  improvement  using 
its  noise  estimate.  The  noise  estimates  of  both  the  AR  related  techniques  were  rather  weak, 
identifying  noise  poles  to  only  within  a  quadrant  of  the  complex  plane.  Thus,  it  is  likely  the  2D 
IQML  estimates  were  already  accurate  enough  that  no  benefit  was  gained  from  the  weak  noise 
estimates.  The  performance  of  the  techniques  in  separable  noise  was  slightly  better  with  respect 
to  the  CRB.  Of  the  2-D  unknown  noise  techniques  2D  UCIQML  improves  over  the  performance  of 
2D  IQML  in  colored  noise,  and  2D  UCMODE  performs  better  then  2D  ARMODE.  The  2D  MODE 
and  related  colored  noise  techniques  fall  short  of  the  full  2D  CRB,  as  they  are  expected  to  do,  since 
they  are  1-D  by  1-D  techniques  [52]. 
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Mean  Squared  Error  of  arg(gamma)  mode  2 


Figure  29.  Estimation  error  for  one  of  two  exponentials  in  exponential  noise  (Stochastic  ML). 

Techniques:  2-D  method  of  direction  estimation  (mode2d),  known  colored  noise 
mode2d  (cmode2d),  unknown  colored  noise  mode2d  (ucmode2d),  and  exponential 
noise  ucmode2d  (armode2d) 


Mean  Squared  Error  of  arg(gamma)  mode  2 


Figure  30.  Estimation  error  for  one  of  two  exponentials  in  separable  noise  (Stochastic  ML).  Tech¬ 
niques:  2-D  method  of  direction  estimation  (mode2d),  known  colored  noise  mode2d 
(cmode2d),  unknown  colored  noise  mode2d  (ucmode2d),  and  exponential  noise  uc- 
mode2d  (armode2d) 
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VII.  Application  to  Synthetic  Aperture  Radar 

In  this  chapter,  the  1-D  and  2-D  colored  noise  techniques  developed  in  the  previous  chapters 
are  applied  to  simulated  and  radar  chamber  SAR  data.  As  a  measure  of  merit,  the  residual  error 
between  the  estimated  model  and  the  data  is  used  to  determine  how  well  the  model  fits  the  data. 
This  metric  is  used  to  determine  if  the  colored  noise  algorithms  provide  better  estimates  of  expo¬ 
nential  content  of  the  data  than  their  white  noise  counterparts.  In  addition,  the  scatterer  locations 
identified  by  the  2-D  techniques  are  examined  for  suitability  as  pattern  recognition  features. 

7.1  1-D  SAR  Data 

To  test  the  algorithms  on  1-D  radar  data,  two  down-range  profiles  (cross-range  data  was  less 
well  modeled  as  damped  exponentials)  were  selected  from  the  SAR  image  of  a  tank  produced  by 
XPATCH  software  [1].  The  profiles  are  shown  in  Figures  31  and  34.  A  damped  exponential  model 
was  fit  to  each  profile  with  colored  and  white  noise  algorithms  for  model  orders  p  =  1 ...  10.  A  noise 
model  order  of  g  =  p,  or  g  =  2  was  used  for  the  unknown  colored  noise  algorithms.  The  normalized 
error  between  the  estimated  model,  y,  and  the  radar  profile  y,  was  calculated  as 

I|y-y|l2 

I|y|l2 

where  ||y||2  =  Vy*y.  The  resulting  modeling  error  is  shown  in  Figures  32  through  36.  The  results 
show  that  for  low  model  order  p,  the  ARUCIQML  and  ARMAIQML  algorithms  fit  the  data  better. 
As  the  model  order  increases,  IQML  and  UCIQML  continually  model  more  of  the  exponential 
signal  and  noise,  and  better  fit  the  data.  Use  of  ARUCIQML  and  ARMAIQML  fits  this  radar  data 
better  than  IQML  at  low  model  orders  and  may  provide  mode  accurate  relative  scatterer  locations 
(exponential  frequencies)  for  target  identification  in  this  case.  Results  for  the  stochastic  algorithm 
for  a  noise  model  order  of  g  =  p  are  shown  in  Figures  37  and  38. 

As  was  indicated  in  the  development  of  the  1-D  colored  noise  algorithms,  there  is  an  easy  way 
to  increase  the  accuracy  of  white  noise  algorithms.  Simply  fit  a  larger  model  order  to  the  data  than 
the  expected  number  of  signals.  Then,  pick  the  highest  energy  part  of  the  estimate  as  signal.  This 
is  called  overmodeling.  It  is  also  the  initial  step  in  the  AR  and  ARMA  colored  noise  algorithms. 
Including  the  overmodeling  technique,  IQMLO,  the  resulting  modeling  error  is  shown  in  Figures 
39  and  40.  A  noise  model  order  of  g  =  p  was  used  for  the  unknown  colored  noise  algorithms. 
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Down-Range  (feet) 


Figure  31.  Down-range  radar  profile  1 


Error  in  Modeling  Radar  Profile  #1 .  q=p 


Figure  32.  Modeling  error  using  p  noise  poles  (Deterministic  ML,  Profile  1).  Techniques:  iterative 
quadratic  maximum  likelihood  (iqml),  unknown  colored  noise  iqml  (uciqml),  autore¬ 
gressive  uciqml  (aruciqml),  and  autoregressive  moving  average  uciqml  (armaiqml) 
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Error  in  Modeling  Radar  Profile  #1 ,  q=2 


Figure  33.  Modeling  error  using  two  noise  poles  (Deterministic  ML,  Profile  1).  Techniques:  iter¬ 
ative  quadratic  maximum  likelihood  (iqml),  unknown  colored  noise  iqml  (uciqml),  au¬ 
toregressive  uciqml  (aruciqml),  and  autoregressive  moving  average  uciqml  (armaiqml) 


Down-Range  Profile  of  a  Tank,  #2 


Figure  34.  Down-range  radar  profile  2 
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Error  in  Modeling  Radar  Profile  #2,  q=p 


Figure  35.  Modeling  error  using  p  noise  poles  (Deterministic  ML,  Profile  2).  Techniques:  iterative 
quadratic  maximum  likelihood  (iqml),  unknown  colored  noise  iqml  (uciqml),  autore¬ 
gressive  uciqml  (aruciqml),  and  autoregressive  moving  average  uciqml  (armaiqml) 


Error  in  Modeling  Radar  Profile  #2,  q=2 


Figure  36.  Modeling  error  using  two  noise  poles  (Deterministic  ML,  Profile  2).  Techniques:  iter¬ 
ative  quadratic  maximum  likelihood  (iqml),  unknown  colored  noise  iqml  (uciqml),  au¬ 
toregressive  uciqml  (aruciqml),  and  autoregressive  moving  average  uciqml  (armaiqml) 
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Error  in  Modeling  Radar  Profile  #1 ,  q=p 
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Model  Order 

Figure  37.  Modeling  error  using  p  noise  poles(Stochastic  ML,  Profile  1).  Techniques:  method  of 
direction  estimation  (mode),  unknown  colored  noise  mode  (ucmode),  autoregressive 
ucmode  (arucmode),  and  autoregressive  moving  average  ucmode  (armamode) 


Error  in  Modeling  Radar  Profile  #2,  q=p 


Model  Order 


Figure  38.  Modeling  error  using  p  noise  poles(Stochastic  ML,  Profile  2).  Techniques:  method  of 
direction  estimation  (mode),  unknown  colored  noise  mode  (ucmode),  autoregressive 
ucmode  (arucmode),  and  autoregressive  moving  average  ucmode  (armamode) 
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Error  in  Modei^g  Radar  Profile  #1 ,  q=p 
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Figure  39.  Modeling  error  using  p  noise  poles  (Overmodeling,  Profile  1).  Techniques:  iterative 
quadratic  maximum  likelihood  (iqml),  iqml  with  p  extra  noise  poles  (iqmlo),  unknown 
colored  noise  iqml  (uciqml),  autoregressive  uciqml  (aruciqml),  and  autoregressive  mov¬ 
ing  average  uciqml  (armaiqml) 


When  overmodeling  is  compared  to  the  colored  noise  techniques  the  performance  difference  is  less 
significant.  In  this  sense,  the  colored  noise  algorithms  may  simply  be  correcting  for  a  bad  choice 
of  model  order.  From  an  algorithmic  point  of  view,  the  difference  between  overmodeling  and  the 
colored  noise  approaches  is  clear.  Both  approaches  estimate  the  noise.  Overmodeling  discards  the 
noise  part  of  the  estimate.  Conversely,  the  colored  noise  techniques  use  the  noise  estimates  to 
improve  the  signal  estimates. 


7.2  2-D  SAR  Data 

For  2-D  data,  bounds  on  the  expected  performance  of  any  technique  that  produces  good 
estimates  are  described  below.  These  bounds  assume  that  the  data  are  well  modeled  as  a  finite  sum 
of  damped  exponentials.  For  the  Synthetic  Aperture  Radar  (SAR)  problem,  the  Mi  x  M2  focused 
image,  X,  and  interpolated  phase  history,  W,  are  represented  as  a  discrete  Fourier  transform  (DFT) 
pair  by 

X  =  DmiWD]^^,  (208) 
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Figure  40.  Modeling  error  using  p  noise  poles  (Overmodeling,  Profile  2).  Techniques:  iterative 
quadratic  maximum  likelihood  (iqml),  iqml  with  p  extra  noise  poles  (iqmlo),  unknown 
colored  noise  iqml  (uciqml),  autoregressive  uciqml  (aruciqml),  and  autoregressive  mov¬ 
ing  average  uciqml  (armaiqml) 

where  Dm  is  the  normalized  DFT  matrix  with  element  0  <  k,l  <  M  —  1.  A 

measure  of  the  error  between  the  phase  histories,  W,  and  the  estimated  model,  W ,  is  given  by 

e  =  \\E\\F  =  \\W-W\\F.  (209) 

To  bound  the  performance  of  parameter  estimation  schemes  one  might  consider  low  rank  models 
other  than  the  damped  exponential  model.  A  general  low  rank  decomposition  of  the  Mi  x  Mg 
matrix,  W,  is 

W  =  Wp  +  E,  (210) 

where  Wp  is  rank-p  and  E  is  a  residual  matrix.  The  best  approximation  of  order  p  is  found  by 
using  the  SVD  of  IT  :  W  =  UHV*.  Here  the  approximated  matrix,  W,  can  be  assembled  from 
rank-p  versions  of  the  SVD  component  matrices  as 

W^Uil:  p,  :)E(1 :  p,  1  :  p)V(l  :  p,  :)*,  (211) 
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where  E  contains  the  p  largest  singular  values,  the  p  columns  of  U  are  the  associated  left  singular 
vectors,  and  the  p  columns  of  V*  are  the  associated  right  singular  vectors.  The  rank-p  SVD 
approximation  provides  a  lower  bound  on  the  error  of  all  other  rank-p  approximations  in  the  form 
W  =  GSH^ .  The  error  attained  by  the  rank-p  SVD  approximation  is  given  by 


esvD  =  ||£'5Vd||f  = 


1 


min  (Mi,Af2) 


*=p+l 


(212) 


where  ctj  is  the  singular  value  in  S.  Now,  consider  the  rank-p  model  of  2-D  damped  exponentials. 
The  columns  of  G  and  H  are  Vandermonde.  Since  this  is  a  restriction,  this  limits  the  accuracy  of 
our  rank-p  approximation.  The  error  incurred  by  this  assumption  is  less  than  that  of  Equation  212. 
Thus  edamp  >  csvD-  Now,  consider  a  further  restriction  on  the  structure  of  G  and  H  where 
the  exponential  frequencies  A,,  and  7,-  can  take  on  only  discrete  values,  specifically,  the  Discrete 
Fourier  Transform  (DFT)  frequencies.  A,-  6  and  7,-  G  Observe  from 

Equations  208  and  209  that  the  minimum  is  attained  by  selecting  the  spatial  frequencies  of  the  p 
values  of  X  with  largest  moduli.  This  DFT  frequency  assumption  represents  a  restriction  to  the 
damped  exponential  model.  The  corresponding  errors  then  satisfy  eopT  ^  ^damp  >  ^SVD- 


7.2.1  Model  Fitting.  The  2-D  exponential  estimation  algorithms  were  tested  with  two 
sources  of  SAR  data  whose  images  are  shown  in  Figures  6,  7,  41,  and  42.  The  first  source  is 
the  XPATCH  radar  cross  section  (RCS)  modeling  tool.  XPATCH  computes  polar-coordinate  RCS 
measurements  from  a  Computer  Aided  Design  (CAD)  model  of  a  target.  The  data  selected  is 
of  a  generic  tank  seen  from  above  at  a  steep  look-angle.  The  data  was  processed  with  the  SAR 
simulation  software  provided  with  XPATCH,  which  produced  images  interpolated  to  a  rectangular 
grid.  Since  the  entire  extent  of  the  data  is  interpolated,  the  phase  history  contains  regions  of  zero 
measurement.  The  windowing  effect  of  including  these  zero  measurements  in  the  data  degrades 
exponential  estimation.  Thus,  a  non-zero  square  segment  of  the  phase  history  data  was  selected. 

The  second  source  of  SAR  data  is  radar  test  chamber  data  of  a  scale  model  of  a  C-29  aircraft. 
Algorithms  were  tested  on  this  unfocused  data  and  on  the  data  focused  by  two  of  the  methods 
discussed  in  Chapter  4.  Due  to  computational  considerations,  only  scattering  center  focus  could  be 
applied  to  the  full  128  x  128  data  size.  Also,  due  to  computational  considerations,  the  techniques 
based  on  2D  IQML  could  only  be  run  on  a  16  x  16  data  size  (a  comparable  computational  load  to 
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Figure  41.  C-29  focused  with  method  from  [42]  using  500  x  500  samples 


Figure  43.  XPATCH  tank,  16  x  16 


running  the  techniques  based  on  2D  MODE  on  the  128  x  128  data  size).  Because  of  the  better  than 
Fourier  resolution,  the  scatterer  locations  of  the  2-D  techniques  are  compared  to  (overlaid  on)  an 
image  of  the  full  128  x  128  data  size.  Two  method  were  used  to  create  the  focused  16  x  16  data.  In 
one  case,  the  128  x  128  scattering  center  focus  image  was  decimated  by  8,  and,  in  the  other  case, 
the  inverse  distance  method  was  used  on  a  16  x  16  sample  of  the  chamber  data.  The  16  x  16  data 
samples  for  the  XPATCH  tank  and  C-29  aircraft  are  shown  in  Figures  43  through  45.  An  example 
of  the  Fourier  resolution  for  16  x  16  data  samples  (zero-padded  in  the  phase  history  to  128  x  128 
samples)  is  shown  in  Figure  46  for  the  XPATCH  tank. 

For  this  set  of  SAR  data,  the  fit  of  the  low  rank  approximations  discussed  in  the  previous 
section  was  determined.  The  normalized  error  between  an  estimated  model  and  the  phase  history 
was  calculated  as 

iniF 

The  error  bounds  given  by  the  SVD  and  the  DFT  were  established  for  this  data.  These  represent 
lower  and  upper  bounds,  respectively,  between  which  any  reasonable  approximation  method  should 
perform.  The  2-D  white  and  colored  noise  techniques  developed  in  the  previous  chapter  were  tested 
with  this  data.  Results  for  the  images  are  plotted  for  the  deterministic  techniques  in  Figures  47 
through  50,  for  the  stochastic  techniques  in  Figures  51  through  53  and  56  through  58  ,  and  for  other 
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Figure  44.  C-29  aircraft,  scatter  center  focus 
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Figure  45.  C-29  aircraft,  inverse  distance  focus 
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Figure  46.  XPATCH  tank,  16  x  16  zero-padded  to  128  x  128 

2-D  techniques  in  Figure  54,  55,  59  and  60.  These  results  compare  how  well  the  exponential  content 
chosen  by  each  technique  fits  the  data.  A  better  fit  will  result  in  less  error.  The  actual  pole  locations 
are  shown  later  for  some  sample  model  orders  to  understand  the  suitability  of  the  techniques  as 
region  of  interest  operators,  or  features  for  model-based  pattern  recognition.  Locations  of  the 
scatterers  in  the  images  are  plotted  for  the  deterministic  techniques  in  Figures  61  through  64,  for 
the  stochastic  techniques  in  Figures  65  through  67  and  70  through  72,  and  for  other  2-D  techniques 
in  Figures  68,  69,  73  and  74. 

The  pairing  necessary  in  the  1-D  by  1-D  techniques  may  be  accomplished  in  several  ways, 
depending  on  an  assumption  of  exponentials  with  distinct  or  non-distinct  frequencies.  The  1-D  by 
1-D  techniques  developed  in  Chapter  6  use  distinct  exponentials  for  an  initial  estimate,  then  allow 
the  final  result  to  have  non-distinct  exponentials.  This  combination  results  in  the  best  performance 
in  terms  of  the  CRB  and  model  fitting  error.  A  drawback  of  this  combination,  the  tendency  of  the 
1-D  by  1-D  methods  to  pick  false  scatterers  in  the  same  column  or  row  as  true  scatterers  (Figure  72), 
can  be  mitigated  by  other  combinations.  Other  choices  of  distinct  and  non-distinct  scatters  can  lead 
to  the  identification  of  additional  scattering  centers  as  shown  in  Figure  75  and  76.  One  drawback 
is  these  cases  is  in  the  fully  distinct  pairing  method;  erroneous  combinations  of  exponentials  will 
result  from  the  forced  pairing.  One  possible  remedy  of  the  pairing  problem  is  to  overmodel  the 
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exponentials  as  in  the  2D  ARMODE  technique  then  choose  the  signal  exponentials  with  some  full 
2-D  matching  scheme.  A  similar  approach  is  taken  in  the  2D  Prony  matching  technique  [59]. 

The  matching  scheme  of  2D  Prony  estimates  the  amplitudes  of  each  row  (column)  based  on 
the  exponential  frequency  estimates  for  the  columns  (rows).  The  amplitudes  are  then  used  as  data 
to  attain  p  row  exponential  frequency  estimates  to  pair  with  the  column  (row)  frequency.  The 
pairings  of  exponentials  closest  to  the  pairings  generated  when  the  roles  of  row  and  colnmn  are 
reversed  (in  parentheses)  are  matched.  A  matching  technique  for  the  1-D  by  1-D  algorithms  along 
the  lines  of  the  2D  Prony  technique  should  at  least  duplicate  the  good  performance  of  this  technique 
in  selecting  many  individual  scatterers  and  not  simply  overmodeling  high-energy  scatterers  (Figures 
73  and  74).  This  is  indeed  the  case  as  shown  in  Figures  77  and  78.  The  fit  of  the  data  using  the 
1-D  by  1-D  stochastic  ML  techniques  is  also  improved  by  matching  in  this  manner  (Figures  79  and 
80).  The  scattering  center  locations  identified  by  the  2-D  deterministic  ML  techniques  are  also 
improved  by  performing  the  estimation  twice  and  matching  the  results  as  in  2D  Prony  (Figures  81 
through  84). 

7.3  Discussion 

Although  the  data  is  limited  some  observations  can  be  made  from  the  results.  The  error 
in  fitting  the  data  for  the  2-D  techniques  do  not  in  general  indicate  better  performance  at  low 
model  as  the  1-D  results  do.  This  is  probably  due  to  an  inability  to  model  the  scatterers  well  in 
both  dimensions  with  damped  exponentials.  The  cross-range  dimension  created  by  the  synthetic 
aperture  is  not  the  smooth  damped  exponential  data  seen  in  the  down-range  data.  Overall,  the 
colored  noise  techniques  show  no  significant  improvement  in  fitting  the  data  over  the  white  noise 
techniques. 

In  the  area  of  focusing,  the  fit  error  of  the  focused  C-29  data  (Figure  58)  vs.  the  unfocused 
C-29  data  (Figure  57)  and  the  scatterers  identified  (Figures  72  vs.  71)  indicates  the  improvement 
gained  from  focusing.  Additionally,  the  inverse  distance  focusing  (Figures  49  and  52)  does  work 
as  well  with  this  data  as  subsampling  the  scattering  center  focus  data  (Figure  50  and  53)  does. 
Scattering  center  focus  is  also  especially  beneficial  for  the  1-D  by  1-D  techniques  since  it  is  the  only 
computationally  tractable  focusing  method  for  the  full  128  x  128  data. 
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Comparing  the  fit  of  different  techniques  for  the  same  size  data  set  (Figures  47,  51  and  54  and 
Figures  50,  53  and  55  )  and  the  scatterer  locations  identified  (Figures  61,  65  and  68  and  Figures 
64,  67  and  69  )  shows  that  the  full  2-D  techniques  based  on  2D  IQML  fit  the  data  the  best,  the  1-D 
by  1-D  techniques  based  on  2D  MODE  performed  next  best,  while  the  remaining  2-D  techniques 
did  not  fit  the  data  as  well. 

The  high  resolution  of  these  techniques  is  evident  in  the  results.  Notice  that  for  most  tech¬ 
niques  the  major  scatterers  in  the  128  x  128  image  are  accurately  located  despite  use  of  only 
16  X  16  data  samples.  For  the  1-D  by  1-D  techniques,  however,  fewer  false  scattering  locations  were 
identified  with  the  128  x  128  data. 

Consistent  identification  of  the  same  scatterers  is  important  for  target  recognition.  The  2-D 
techniques  based  on  2D  IQML  consistently  identify  the  same  scattering  centers  in  one  image  at 
different  model  orders  (Figure  61).  However,  since  2D  IQML  requires  exponentials  with  distinct 
frequencies  in  one  dimension,  not  all  the  same  scatterers  are  identified  when  the  image  is  rotated 
90  degrees  (Figure  62)  and  a  different  fit  of  the  model  and  data  is  observed  (Figures  47  and  48). 
Scatterer  locations  are  not  distinct  in  both  dimensions,  thus,  they  are  not  correctly  located  by  these 
techniques  in  all  cases.  For  this  SAR  data  then,  the  use  of  2D  Prony  matching  with  the  exponential 
estimation  techniques  gives  better  scatterer  locations  (Figures  85  and  86). 

The  colored  noise  techniques  with  2D  Prony  matching  produce  a  better  fit  to  the  data  than 
the  white  noise  techniques  at  many  model  orders.  In  addition,  the  colored  noise  techniques  do 
identify  different  scattering  centers  than  the  white  noise  techniques.  The  1-D  by  1-D  techniques, 
as  indicated  earlier,  tend  to  identify  scatterers  on  a  line.  Additionally,  on  the  128  x  128  data  the 
colored  techniques  tended  to  pick  out  new  scattering  centers  instead  of  overmodeling  high  energy 
ones.  The  scattering  centers  identified,  however,  are  dependent  on  assumptions  of  distinct  and 
non-distinct  scattering  centers,  and  many  false  scattering  centers  are  identified  as  discussed  in  the 
previous  section.  This  is  then  the  reverse  of  the  problem  of  2D  IQML  where  too  few  scattering 
centers  are  identified  because  of  the  technique’s  requirement  for  distinct  frequencies  in  one  of  the 
dimensions.  The  solution  to  this  problem  is  provided  by  the  matching  technique  in  the  2D  Prony 
algorithm.  With  this  matching  technique  most  scattering  centers  are  identified  and  false  scattering 
centers  off  the  target  are  eliminated. 
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Figure  47.  Representation  error  vs.  model  order  (Deterministic  ML,  16  x  16  tank  image).  Tech¬ 
niques;  2-D  iterative  quadratic  maximum  likelihood  (2D  IQML),  unknown  colored 
noise  iqml2d  (2D  UCIQML),  and  exponential  noise  uciqml2d  (2D  ARIQML),  Singu¬ 
lar  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT).  2D  IQML  related 
techniques  are  limited  to  estimating  8  total  exponentials  and  noise  poles  (4  exponen¬ 
tials  and  4  noise  poles  in  2D  ARIQML  here) 


Tank 


Figure  48.  Representation  error  vs.  model  order  (Deterministic  ML,  16  x  16  rotated  tank  image). 

Techniques:  2-D  iterative  quadratic  maximum  likelihood  (2D  IQML),  unknown  colored 
noise  iqml2d  (2D  UCIQML),  and  exponential  noise  uciqml2d  (2D  ARIQML),  Singular 
Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 
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Figure  49.  Representation  error  vs.  model  order  (Deterministic  ML,  16  x  16  C-29  inverse  distance 
focus  image).  Techniques;  2-D  iterative  quadratic  maximum  likelihood  (2D  IQML), 
unknown  colored  noise  iqml2d  (2D  UCIQML),  and  exponential  noise  uciqml2d  (2D 
ARIQML),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 


c-29 


Figure  50.  Representation  error  vs.  model  order  (Deterministic  ML,  16x16  C-29  scattering  center 
focus  image).  Techniuqes:  2-D  iterative  quadratic  maximum  likelihood  (2D  IQML), 
unknown  colored  noise  iqml2d  (2D  UCIQML),  and  exponential  noise  uciqml2d  (2D 
ARIQML),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 
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Figure  51.  Representation  error  vs.  model  order  (Stochastic  ML,  16  x  16  tank  image).  Techniques; 

2-D  method  of  direction  estimation  (2D  MODE),  unknown  colored  noise  mode2d  (2D 
UCMODE),  and  exponential  noise  ucmode2d  (2D  ARMODE),  Singular  Value  Decom¬ 
position  (SVD),  Discrete  Fourier  Transform  (DFT) 


Tank 


Figure  52.  Representation  error  vs.  model  order  (Stochastic  ML,  16  x  16  C-29  inverse  distance 
focus  image).  Techniques:  2-D  method  of  direction  estimation  (2D  MODE),  un¬ 
known  colored  noise  mode2d  (2D  UCMODE),  and  exponential  noise  ucmode2d  (2D 
ARMODE),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 
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Figure  53.  Representation  error  vs.  model  order  (16  x  16  C-29  scattering  center  focus  image). 

Techniques:  2-D  method  of  direction  estimation  (2D  MODE),  unknown  colored  noise 
mode2d  (2D  UCMODE),  and  exponential  noise  ucmode2d  (2D  ARMODE),  Singular 
Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 


Figure  54.  Representation  error  vs.  model  order  (Other  techniques,  16  x  16  tank  image).  Tech¬ 
niques:  2-D  Prony  (2D  PRONY),  state  space  (STATESP),  and  matrix  enhance¬ 
ment  matrix  pencil  (MEMP)  methods.  Singular  Value  Decomposition  (SVD),  Discrete 
Fourier  Transform  (DFT) 
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Figure  55.  Representation  error  vs.  model  order  (Other  techniques,  16  x  16  C-29  scattering  center 
focused  image).  Techniques:  2-D  Prony  (2D  PRONY),  state  space  (STATESP),  and 
matrix  enhancement  matrix  pencil  (MEMP)  methods,  Singular  Value  Decomposition 
(SVD),  Discrete  Fourier  Transform  (DFT) 


Tank 


Figure  56.  Representation  error  vs.  model  order  (Stochastic  ML,  128  x  128  tank  image).  Tech¬ 
niques:  2-D  method  of  direction  estimation  (2D  MODE),  unknown  colored  noise 
mode2d  (2D  UCMODE),  and  exponential  noise  ucmode2d  (2D  ARMODE),  Singu¬ 
lar  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 
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Figure  57.  Representation  error  vs.  model  order  (Stochastic  ML,  128x128  C-29  unfocused  image). 

Techniques:  2-D  method  of  direction  estimation  (2D  MODE),  unknown  colored  noise 
mode2d  (2D  UCMODE),  and  exponential  noise  ucmode2d  (2D  ARMODE),  Singular 
Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 


c-29 


Figure  58.  Representation  error  vs.  model  order  (Stochastic  ML,  128  x  128  C-29  scattering  cen¬ 
ter  focused  image).  Techniques:  2-D  method  of  direction  estimation  (2D  MODE), 
unknown  colored  noise  mode2d  (2D  UCMODE),  and  exponential  noise  ucmode2d  (2D 
ARMODE),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 
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Figure  59.  Representation  error  vs.  model  order  (Other  techniques,  64  x  64  tank  image).  Tech¬ 
niques:  2-D  Prony  (2D  PRONY),  state  space  (STATESP),  and  matrix  enhance¬ 
ment  matrix  pencil  (MEMP)  methods,  Singular  Value  Decomposition  (SVD),  Discrete 
Fourier  Transform  (DFT) 


C-29 


Figure  60.  Representation  error  vs.  model  order  (Other  techniques,  64  X  64  C-29  scattering  center 
focused  image).  Techniques:  2-D  Prony  (2D  PRONY),  state  space  (STATESP),  and 
matrix  enhancement  matrix  pencil  (MEMP)  methods.  Singular  Value  Decomposition 
(SVD),  Discrete  Fourier  Transform  (DFT) 
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Model  Order  =  7 


Figure  61.  Deterministic  ML  scatterers  (16x16  tank  image). 

o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 


Model  Order  =  8 


Figure  62.  Deterministic  ML  scatterers  (16  x  16  rotated  tank  image), 
o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 


Model  Order  =  8 


Figure  63.  Deterministic  ML  scatterers  (16  x  16  C-29  inverse  distance  focused  image), 
o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 
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Figure  64.  Deterministic  ML  scatterers  (16  x  16  C-29  scattering  center  focus  image), 
o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 
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Figure  65.  Stochastic  ML  scatterers  (16  x  16  tank  image). 

o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Figure  66.  Stochastic  ML  scatterers  (16  x  16  C-29  inverse  distance  focused  image), 
o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Figure  67.  Stochastic  ML  scatterers  (16  x  16  C-29  scattering  center  focused  image), 
o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Model  Order  =  7 
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Figure  69.  Scattering  centers  of  other  2-D  techniques  (16  x  16  C-29  scattering  center  focused 
image) . 

o  =  2D  PRONY  +  =  STATE  SPACE  x  =  MEMP 
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Figure  71.  Stochstic  ML  scatterers  (128  x  128  C-29  unfocused  image), 
o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Figure  72.  Stochastic  ML  scatterers  (128  x  128  C-29  scattering  center  focused  image), 
o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Model  Order  =  7 


pixel 

Figure  73.  Scattering  centers  of  other  2-D  techniques  (64  x  64  tank  image), 
o  =  2D  PEONY  +  =  STATE  SPACE  x  =  MEMP 
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Model  Order  =  7 


Figure  74.  Scattering  centers  of  other  2-D  techniques  (64  x  64  C-29  scattering  center  focused 
image). 

o  =  2D  PRONY  +  =  STATE  SPACE  x  =  MEMP 
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Figure  75.  Stochastic  ML  distinct  scatterers  (128  x  128  C-29  scattering  center  focused  image), 
o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Figure  76.  Stochastic  ML  non-distinct  scatterers  (128 x  128  C-29  scattering  center  focused  image), 
o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Figure  78.  Stochastic  ML  scatterers,  2D  Prony  matching  (64  x  64  C-29  scattering  center  focused 
image). 

o  =  2D  MODE  +  =  2D  UCMODE  x  =  2D  ARMODE 
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Tank 


Figure  79.  Representation  error  vs.  model  order  (Stochastic  ML,  64  x  64  tank  image),  2D  Prony 
matching.  Techniques:  2-D  method  of  direction  estimation  (2D  MODE),  unknown  col¬ 
ored  noise  mode2d  (2D  UCMODE),  and  exponential  noise  ucmode2d  (2D  ARMODE), 
Singular  Value  Decomposition  (SVD),  Discrete  Fourier  Transform  (DFT) 


C-29 


Figure  80.  Representation  error  vs.  model  order  (Stochastic  ML,  64  x  64  C-29  scattering  center 
focused  image),  2D  Prony  matching.  Techniques:  2-D  method  of  direction  estimation 
(2D  MODE),  unknown  colored  noise  mode2d  (2D  UCMODE),  and  exponential  noise 
ucmode2d  (2D  ARMODE),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier 
Transform  (DFT) 
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Figure  81.  Deterministic  ML  scatterers,  2D  Prony  matching  (16  x  16  Tank  image), 
o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 


120 


129 


Model  Order  =  7 


pixel 


Figure  82.  Deterministic  ML  scatterers,  2D  Prony  matching  (16  x  16  C-29  scattering  center  fo¬ 
cused  image). 

o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 
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Figure  83.  Representation  error  vs.  model  order  (Deterministic  ML,  16  x  16  tank  image), 
2D  Prony  matching.  Techniques:  2-D  iterative  quadratic  maximum  likelihood  (2D 
IQML),  unknown  colored  noise  iqml2d  (2D  UCIQML),  and  exponential  noise  uciqml2d 
(2D  ARIQML),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier  Transform 
(DFT) 


C-29 


Figure  84.  Representation  error  vs.  model  order  (Deterministic  ML,  16  x  16  C-29  scattering  center 
focused  image),  2D  Prony  matching.  Techniques:  2-D  iterative  quadratic  maximum 
likelihood  (2D  IQML),  unknown  colored  noise  iqml2d  (2D  UCIQML),  and  exponential 
noise  uciqml2d  (2D  ARIQML),  Singular  Value  Decomposition  (SVD),  Discrete  Fourier 
Transform  (DFT) 
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Deterministic  ML  scatterers,  2D  Prony  matching  (16  x  16  C-29  scattering  center  fo 
cused  image). 

o  =  2D  IQML  +  =  2D  UCIQML  x  =  2D  ARIQML 


VIII.  Conclusions 


Maximum  likelihood  techniques  were  developed  to  estimate  exponentials  in  colored  noise. 
Methods  were  developed  for  one  and  two  dimensions,  for  stationary  colored  noise,  and  colored 
noise  modeled  with  a  rational  polynomial.  For  all  these  cases,  separate  methods  were  developed 
for  a  deterministic  signal  assumption  and  for  a  stochastic  signal  assumption. 

All  these  techniques  perform  well  when  the  noise  model  matches  the  actual  noise.  In  these 
cases,  accuracy  in  estimating  exponential  frequencies  was  improved  over  techniques  assuming  white 
noise.  All  techniques  attained  the  Cramer-Rao  estimation  bound  (CRB)  either  for  short  data 
records  or  asymptotically  with  good  performance  for  short  data  records.  These  techniques  also 
performed  well  when  the  noise  was  white  and,  in  some  cases,  improved  low  SNR  performance  by 
mitigating  the  colored  appearance  of  a  finite  sample  of  white  noise.  In  two  dimensions,  the  increased 
degrees  of  freedom  of  2-D  data  makes  it  difficult  to  accurately  estimate  the  noise.  Thus,  only  slight 
improvement  over  the  white  noise  techniques  was  observed  and  no  technique  attained  the  CRB  for 
short  data  records. 

The  colored  noise  techniques  were  applied  to  Synthetic  Aperture  Radar  data  which,  when 
correctly  processed,  can  be  modeled  as  the  sum  of  damped  exponentials  for  targets  consisting  of 
point  scatterers.  In  one  dimension,  colored  noise  techniques  corrected  for  a  poor  choice  of  model 
order  in  much  the  same  way  as  modeling  the  data  with  an  increased  model  order  and  discarding 
the  low  energy  part  (overmodeling).  In  fact,  the  I-D  colored  noise  estimation  techniques  can  be 
viewed  as  overmodeling  then  using  the  overmodeled  part  as  a  noise  estimate  and  to  improve  the 
estimate  of  the  signal.  In  two  dimensions,  better  fit  of  the  data  at  low  model  orders  was  not 
observed.  The  2-D  techniques  based  on  2D  IQML  identified  a  limited  number  of  scatterer  in  one 
dimension  due  to  a  limitation  on  the  algorithm  to  distinct  frequencies  in  one  dimension.  The  1-D 
by  1-D  techniques  based  on  2D  MODE  tended  to  identify  scatterers  along  lines  in  either  dimension. 
The  use  of  2-D  Prony  matching  improved  scattering  center  locations  in  all  techniques.  The  2-D 
techniques  based  on  2D  IQML  and  on  2D  MODE  when  paired  with  the  matching  techniques  of  2D 
Prony  consistently  identified  many  individual  scatterers  and  had  few  false  scatterers.  In  these  cases, 
the  unknown  colored  noise  techniques  identified  different  scatterers  from  the  white  noise  techniques 
and  fit  the  data  better  as  well. 
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The  2-D  techniques  based  on  the  colored  noise  model  more  accurately  model  SAR  data 
than  existing  2-D  white  noise  techniques.  The  2-D  colored  noise  techniques  developed  in  this 
thesis  should  significantly  aid  in  the  development  of  model-based  pattern  recognition  algorithms  to 
identify  targets  in  SAR  images.  All  of  the  2-D  techniques  when  combined  with  scattering  center 
focus  and  2D  Prony  matching  provide  robust  tools  for  identifying  scattering  centers  on  an  object 
in  a  SAR  image.  These  techniques  will  aid  in  understanding  the  consistency  with  which  scatterers 
can  be  identified  across  different  SAR  images.  Of  the  techniques  the  2D  IQML  and  2D  UCIQML 
techniques  tied  with  2D  Prony  matching  appear  to  give  the  best  estimates  of  the  locations  of  true 
2-D  scatterers  while  avoiding  overmodeling  of  these  scatterers. 

8.1  Contributions 

List  of  the  contributions  in  this  dissertation. 

•  A  concise  mathematical  model  for  damped  exponentials  on  an  irregularly  sampled  grid. 
(Chapter  4) 

•  A  computationally  efficient  (real-time)  method  for  interpolating  or  focusing  SAR  images 
containing  damped  exponentials.  (Chapter  4) 

•  Closed-form  solution  of  the  1-D  maximum  likelihood  problem  of  estimating  deterministic 
exponentials  in  unknown  colored  noise.  (Chapter  5) 

•  Reformulation  of  the  1-D  colored  noise  stochastic  maximum  likelihood  problem  to  avoid  a 
global  search  for  the  solution.  (Chapter  5) 

•  Several  1-D  deterministic  and  stochastic  maximum  likelihood  algorithms  for  estimating  1-D 
exponentials  in  colored  noise.  (Chapter  5) 

•  Extension  of  deterministic  maximum  likelihood  estimation  of  2-D  exponentials  to  unknown 
colored  noise.  (Chapter  6) 

•  Extension  of  stochastic  maximum  likelihood  estimation  of  2-D  exponentials  on  a  1-D  by  1-D 
grid  to  unknown  colored  noise.  (Chapter  6) 

•  A  spectral  model  to  produce  unique  2-D  spectra  (noise,  or  an  impulse  driving  a  2-D  damped 
exponential  filter).  (Chapter  6) 
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•  Several  2-D  deterministic  and  stochastic  maximum  likelihood  algorithms  for  estimating  2-D 
exponentials  in  colored  noise  (Chapter  6) 

•  Unique  properties  of  lower  (upper)  triangular  Toeplitz  matrices,  commutivity  and  a  lower 
(upper)  triangular  Toeplitz  inverse.  (Appendix  A) 

8.2  Future  Work 

There  are  several  areas  where  the  work  initiated  by  this  dissertation  could  be  expanded  and 
potential  new  insights  gained; 

•  Apply  the  unknown  colored  noise  techniques  developed  in  this  dissertation  to  multiple  data 
instances  and  unit-circle  constrained  frequencies.  In  this  case,  many  other  application  are 
available  for  the  techniques,  such  as,  array  processing  in  the  presence  of  noisy  jammers. 

•  Develop  methods  of  simultaneously  estimating  the  signal  and  the  noise  for  the  unknown  col¬ 
ored  noise  exponential  estimation  methodology  developed  in  this  dissertation.  Such  methods 
will  likely  lead  to  efficient  techniques  with  respect  to  the  CRB  for  exponential  and  noise 
estimation. 

•  Investigate  the  consistency  and  resolution  with  which  scatterers  are  identified  at  different 
look  angles  to  the  target.  Recognition  of  different  targets  requires  consistent  identification  of 
the  scatterer  locations  on  the  targets.  The  large  set  of  XPATCH  (three  SAR  bands,  several 
targets,  several  elevations  and  all  Eispects),  and  C-29  chamber  data  (many  elevations  and 
many  aspects  at  very  fine  frequency  and  angle  resolution)  allows  for  an  extensive  analysis. 

•  Develop  a  model-based  pattern  recognition  system  using  relative  locations  of  scatterers  on  the 
target.  Investigation  of  real  world  effects  such  as  occlusion  of  scattering  centers  will  determine 
the  utility  of  such  a  system. 
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Appendix  A.  Commutivity  of  Operations 


The  approximate  commutivity  of  the  operations  in  Equation  105  is  based  on  the  commutivity 
of  the  convolution  operation.  Similarly,  the  matrix  formulation  of  convolution,  multiplication  of 
lower  (upper)  triangular  Toeplitz  matrices  commutes.  The  commutivity  of  block-diagonal  triangular 
Toeplitz  matrices  is  implied  in  [22].  The  following  theorems  establish  the  commutivity  of  the 
triangular  Toeplitz  matrices,  the  triangular  Toeplitz  structure  of  the  inverse  of  a  triangular  Toeplitz 
matrix,  and  the  commutivity  of  the  operations  in  Equation  105. 

Theorem.  Let  A  ,  B  he  N  x  N  upper  (lower)  triangular  Toeplitz  matrices  and  their  product 
be  C  =  AB,  then,  the  matrices  A  and  B  commute,  and  C  =  BA.  Proof.  The  elements  of  the 
product  of  the  upper  triangular  matrices  A  and  B  are 

J 

Cij  —  i  ^  j.  (214) 

k—t 

Since  the  matrices  A  and  B  are  also  Toeplitz,  they  are  symmetric  about  the  anti-diagonal  and  the 
their  elements  are  related  as 

0'i,k={i...j]  —  ^k={j...i},j'  (215) 


Thus, 

j 

Cij  —  'y  \akj  k  i  ^  i)  (216) 

k~i 

or  C  =  BA.  The  lower  triangular  case  is  simply  seen  by  taking  the  transpose  of  the  upper  triangular 
matrix  products  AB  =  BA. 

Theorem.  Let  T  be  an  upper  (lower)  triangular  Toeplitz  matrix  then  its  inverse  T~^  is  also 
upper  (lower)  triangular  Toeplitz.  Proof.  The  general  form  for  the  inverse  of  a  Toeplitz  martix 
given  in  [45]  is  factored  into  Toeplitz  matrices  as 


(217) 
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where 


and 


ra  = 


Th  = 
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(218) 


(219) 


If  T  is  upper  triangular,  t,-  =  0  i  =  —1 —  N,  and 
simply  the  first  column  of  the  inverse,  normalized  as  a  = 
Equation  217  gives 


since  T  ^ 

[  1  0  . 


is  also  upper  triangular  a  is 


0 


.  Substituting  this  in 


pT-^ 


I  bi  ■■■  bN 
0  : 

:  •••  •••  6i 

0  •••  0  1 


(220) 


Thus  T“^  is  Toeplitz.  The  lower  triangular  case  follows  similarly. 


Theorem.  The  operations  in  Equation  105  commute.  Proof.  Let  Z  be  the  N  x  N  —p  matrix 
consisting  of  the  N  —p  identity  matrix  augmented  with  p  rows  of  zeros  and  Z  be  the  N  x  N  matrix 
consisting  of  Z  augmented  with  p  columns  of  zeros.  Then  A  is  the  N  x  N  lower  triangular  Toeplitz 
matrix  such  A  =  AZ.  The  inverse  of  the  covariance  matrix  R.~^  =  is  approximated 

by  zeroing  the  last  p  rows  of  B  of  with  K  —  ZB  and  resizing  B  with  K  =  Z*  B.  The  neglected  terms 
constitute  an  approximation  with  0{p/N)  error.  The  operations  in  Equation  105  are  approximated 
as 

{A*{BP-^B*)-^A)-^  ~  {A*{KK*)^A)-^  (221) 


For  q  <  p  ,  excepting  zeroed  elements  K  =  ZB  is  lower  triangular  Toeplitz.  Then  excepting  zeroed 
elements,  matrices  (K)'^  and  A  =  ZA  are  lower  triangular  Toeplitz  and  the  same  rank  and  thus 


138 


they  commute. 


{A*{KK*YA)-^  =  {Z*A*Z{KK*)-^ZAZ)-^ 

(222) 

(223) 

=  {Z*  {!<*)+ A*  A{k)-^Z)-^ 

(224) 

=  k{A*A)+k* 

(225) 

since  ZZ  —  Z  and  Z"*"  =  Z* .  The  {A*  A)'^  term  is  also  approximated  with  0{p/N)  error  by  {A*  A)  ^ 
or  Z{A*A)-^Z*. 
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Appendix  B.  Consistent  1-D  Noise  Estimates 

Consistent  estimates  of  stationary,  autoregressive  (AR),  and  autoregressive-moving  average 
(ARMA)  are  attained  through  application  of  the  following  theorems. 

Theorem  (stationary).  A  consistent  estimate  of  the  noise  covariance  matrix  Rww  of  station¬ 
ary  noise  sequence  is  attained  from  the  difference  between  the  observed  data  and  the  maximum 
likelihood  estimate  of  the  signal.  Proof.  Given  the  Toeplitz  structure  of  the  stationary  noise 
covariance  matrix  Rww  and  a  consistent  estimator  of  the  elements  of  the  covariance  matrix  the 
autocorrelation  coefficients,  all  that  remains  to  be  shown  is  that  the  residuals  of  the  maximum 
likelihood  estimate  are  a  consistent  estimate  of  the  noise  sequence.  The  estimate  is 

w  =  y  -  Gs  (226) 

where  Gs  is  the  ML  estimate  of  the  signal.  By  the  triangle  inequality 

||y-Gs-(y-Gs)||  (227) 

||y-y||+||G's-Gs|| 

||Gs-Gs|| 

Since  ||  •  ||  >  0,  squaring  both  sides  and  taking  the  limit  of  the  expected  value  gives 

lim  E{||w  -  w||^}  <  lim  E{||Gs  -  Gs||^}.  (228) 

N—*-oo  AT— *-oo 

And  since  the  ML  estimates  of  the  signal  are  consistent 

^l^E{||w- w|||}  =  0.  (229) 

Theorem  (AR).  Consistent  estimates  of  a  AR(p)  noise  process  are  given  by  the  approxi¬ 
mating  the  maximum  likelihood  minimization  as  shown  in  Appendix  A.  Asymptotically,  the  same 
minimization  gives  consistent  estimates  of  the  AR  noise  process.  Proof.  For  a  fixed  signal  the  ML 
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estimate  of  the  noise  is  given  by 


min-^ln|y*A(^*ii„w^)  =  min In  |y* A(yl*i?ww^)  In  (230) 

■Rww  iV  Ry/uvi  Jy  -A' 

Taking  the  limit  and  introducing  the  reflection  coefficient  decomposition  of  In  |i?ww|  from  [51]  where 
|fc„|  <  1  gives 

1  P 

minlny*A(^*i?ww^)“^^*y  -  hm  —  Y'  n(l  - 

Rww  N— ►oo  iV 

n=l 

=  miny*A(yl*i?ww^)-i^*y.  (231) 

•^w  w 

Theorem  (ARM A).  Consistent  estimates  of  a  ARMA(p)  noise  process  are  also  given  by  the 
approximating  the  maximum  likelihood  minimization  as  shown  in  Appendix  A.  The  AR  and  MA 
parts  of  the  noise  can  be  estimated  separately  when  the  AR  part  of  the  noise  process  contains  most 
of  the  noise  energy  (i.e.  the  noise  poles  are  close  to  the  unit  circle  in  the  complex  plane,  and  the 
noise  poles  are  separated  in  frequency  from  the  moving  average  zeros).  Proof.  The  AR  process 
is  estimated  first.  Then  the  AR  theorem  may  be  applied.  The  stationary  theorem  may  then  be 
applied  to  the  whitened  residuals  of  the  noise  process,  wma  =  where  Rar  =  L^^Lar 

and  (•)'^  indicates  pseudoinverse. 
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Appendix  C.  A  Bound  on  the  MSE  of  Estimates 

r  tT 


A  useful  measure  of  how  well  the  estimates,  0  — 
0,  of  a  model  is  Mean  Square  Error  (MSE), 


^1  02 


,  match  the  parameters. 


MSE  =  (^1  -  0if. 


(232) 


When  the  error  covariance  matrix  is  defined  as 

C^E{i0-0)i0-0f},  (233) 

the  MSE  of  the  parameter  is  the  diagonal  element  of  the  matrix  C  ,  {C)u.  And  when 
the  estimate  is  unbiased,  E{0}  =  0,  the  error  covariance  matrix  is  the  covariance  matrix  of  the 
estimator, 


C  =  E{(^-E{d}  +  E{^}-^)(^-E{^}  +  E{^}-0)^}  (234) 

=  E{(0  -  E{^})(0"  -  E{0})^}+(E{^}  -  0)(E{0}  -  0^  (235) 

=  E{(^-E{0})(^-E{^})^}.  (236) 

In  Maximum  Likelihood  (ML)  estimation,  the  estimates  of  the  parameters  are  given  by  the 
zeros  of  the  score  function  s(d,  x),  the  derivative  of  the  probability  density,  /,  of  the  observed  data, 
X,  with  respect  to  the  parameters  0, 


s(^',x)  =  ^ln/«(x)  =  0.  (237) 

When  the  ML  estimate  exists,  it  is  unbiased  and  a  bound  on  the  ML  MSE  is  the  Cramer-Rao 
Bound  (CRB).  The  following  two  facts  reveal  this  bound. 

1)  The  score  function  is  identically  correlated  with  the  error  [66].  Since  the  estimate  is 
unbiased 

E{0-6»}  =  O  (238) 

or  equivalently 

j  dxf0{x){0  -  0)^  =  0^  (239) 
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Under  regularity  conditions,  the  derivative  with  respect  to  ^  of  this  equation  is  given  by 


j  dx-^fe{x){e  -  J  £^x/e(x)/  =  0^,  (240) 

or  since  jgfe{x.)  =/e(x)^ln/e(x)  and  f  dxfeix)  =1, 

j  dx/9(x)^ln/s(x)(^-  =  I.  (241) 

Which  is  the  first  fact, 

E{s(^,x)(^-^)^}  =  7.  (242) 

2)  A  matrix  inequality  [51].  The  covariance  matrix  of  the  random  vector  y  —  RyzRj^z  is 
positive  semidefinite, 

E{(y  —  7?y^7?^/z)(y  —  Ry^R^J^z)  }  =  Ryy  —  2RyzRj,} R^^  +  Ry^R^^  R^gR^^  Ry^ 

=  ~  ^'yz  —  0-  (243) 

Let  y  =  9  —  9  and  z  =  s(^,  x)  then  Ryz  —  I  from  the  first  fact  and  from  the  second  fact  the 
CRB  is 

C  >  (E{s(^,x)s(^,x)^})-i  =  J-\9),  (244) 

where  J{9)  =E{s(0,  x)s(0,  x)^}  is  the  Fisher  Information  matrix.  Then  the  MSE  of  the  param¬ 
eter  is  bounded  by  the  f**  diagonal  term  of  the  inverse  of  the  FIM. 

(245) 
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Appendix  D.  Cramer-Rao  Bound  for  Deterministic  Exponentials 


This  appendix  describes  the  Cramer-Rao  (CR)  bound  for  damped  exponentials  in  colored 
noise.  The  bounds  are  computed  by  inverting  the  Fisher  information  matrix  (FIM).  The  general 
form  for  the  element  of  the  FIM  for  a  complex  circular  Gaussian  process  [71]  is  described  in 
terms  of  the  mean  Gs  and  covariance  cr^iJvvw  of  the  damped  exponential  in  colored  noise  process 


T  on  _2  „_1  6Gs  _2  0-1  ScB Ryi-w  _2  n-l  Rvrvr  \ 

~  ^ww^)  +  tr(<^  ^ww  ^ww  )> 


where  the  parameters  are  represented  in  the  real  vector 


0  =[  <T^  \r]\'^  arg(7/)^  |s|^  arg(s)^  |A|’’  arg(A)^]^  ,  (247) 

whose  k*^  element  is  Ok-  As  is  the  case  for  undamped  exponentials  [16]  the  estimation  bounds  for 
the  exponentials  and  the  colored  noise  are  decoupled  and  the  FIM  has  the  structure, 

0 

J  — 

0  J(s,A) 

where 

=  (249) 

(250) 

This  decoupling  is  due  to  the  Wold  decomposition  [51]  which  shows  that  any  regular  wide  sense  sta¬ 
tionary  random  process  can  be  separated  into  a  regular  random  process  and  a  purely  deterministic 
process.  Specific  bounds  for  an  AR  noise  process  and  the  77)  block  of  the  FIM  are  developed 
in  [16].  The  structure  of  the  damped  exponential  block  has  also  been  developed  [67],  [10]  for  the 
case  of  white  noise  and  colored  noise,  respectively. 

Compactly  written  the  damped  exponential  bounds  are 

J(s,  A)  =  ^Re{F*R-l,F)  (251) 
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where  the  TV  x  4p  matrix  F  is 


F  =  [GT  iGE  DGT,A-^  jDGT.  ]  (252) 

where 

T  =  diag{  exp(jarg(si))  exp(i arg(s2))  exp(j arg(sp))  },  (253) 

E  =  diag{  Si  S2  •••  Sp  }.  (254) 

A  =  diag{  |Ai|  IA2I  •••  \Xp\  },  (255) 

D  =  diag{  0  1  •••  TV  }•  (256) 

The  CR  bounds  are  then  given  by  E($f.  —  Ok)"^  >  {J~^)kk- 
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Appendix  E.  Method  of  Direction  Estimation  (MODE) 

The  MODE  algorithm  is  described  in  [71]  [69],  and  [82].  MODE  minimizes 

mmiv{A{A*A)-'-A*VWV*)  (257) 

where  the  eigendecomposition  of  the  data  matrix  is  YY*  =  V{D  —  <T^Ipxp)V*  +  and  the 

columns  V  are  the  eigenvectors  of  Y  associated  with  the  p  largest  eigenvalues  (the  elements  of  the 
diagonal  matrix  D),  (E  is  the  average  of  the  remaining  eignevalues  and  the  optimal  weighting  [50] 
is  W  =  {D  —  Ip-^pY .  The  minimization  in  equation  257  is  computed  in  a  similar  manner  to 
IQML  as 

mina*E*Eai.  (258) 

a* 

where  F-[CFi  CFp  7  and  C  =  [A* 


bp+l,n  ■  ■  ■ 

Vp+2,n  •  •  •  V2,„ 


Vm.xi 


Ym —p,n 


and  Vki  is  the  element  of  Essentially  a  data  matrix  similar  to  Y  is  formed  from 

each  eigenvector  in  V  weighted  by  and  these  data  matrices  are  stacked  (for  simultaneous 

minimization). 
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Appendix  F.  A  lower  memory  method  for  implementing  2-D  IQML 

In  the  original  version  of  2D-IQML  the  following  procedure  produce  the  kernel  that  is  used 
to  predict  the  2-D  frequencies: 

T={WW*)^Y 
K  =:  Y*T. 

Since  W  is  ((2Mi  —  1)(M2  -  p)  -|-  Mi  —  l)x  Mi  M2  this  requires  the  SVD  and  pseudoinverse 
of  a  ((2Mi  -  1)(M2  -  p)  -I-  Ml  -  1)  X  ((2Mi  —  1)(M2  —  p)  +  Mi  -  1)  matrix.  The  same  kernel  can 
be  produced  with  less  memory  by  taking  advantage  of  the  smaller  column  width  of  W  and  the  fact 
that  {WW*)^  =  {{W*W)^W*)*{W*W)^W*  which  give  us  the  following  procedure: 

s  =  {w*wYw*y 
K  =  S*S 

We  can  see  that  the  results  are  equivalent  by  expanding  both  equations.  In  the  first  case 
K  =  Y*((UDV*)(UDV*)*)+Y  =  Y*  {U  DV*V  D*U*)'^Y 
=  Y*{UDD*U*)-^Y  =  Y*U{DD*yU*Y 

In  the  second  case 

S  =  {(UDV*)*{UDV*))+iUDVyY  =  {VD*U*UDV*)+VD*U*Y 
=  {VD*DV*)+VD*U*Y  -  V{D*D)+V*VD*U*Y  =  V{D* D)+ D*U*Y . 

And 

K  =  S*S  =  {V{D*D)+D*U*YyV{D*D)+D*U*Y 
=  Y*UD{D*D)+*V*V{D*D)+D*U*Y  =  Y*UD{D*  D)+*  {D*  D)+ D*U*Y 
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where  all  the  inner  matrices  are  diagonal  and  commute.  Thus 


K  =  Y*UiD*D)+*{D*D)+{D*D)U*Y  Y*U{DD*)+U* 


or  the  same  result  as  the  first  case. 
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